Incorrect Series with EllipticTheta under 11.3

Question

The following program gives the correct coefficients under versions 7, 10.2 or 11.0

 Table[SeriesCoefficient[(1/2) x^(-1/4) (EllipticTheta[2, 0, x^1] EllipticTheta[3, 0, x^21] - EllipticTheta[3, 0, x^1] EllipticTheta[2, 0, x^21])/(QPochhammer[x^3] QPochhammer[x^7]), {x, 0, n}], {n, 0, 20}]
 (* {1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 3, 2, 3} *)

But under version 11.3 I got

 (* {1, 0, 1, 1, 0, 1, 3, 1, 2, 5, 1, 3, 9, 3, 7, 13, 6, 9, 21, 10, 18} *)

Why ?

Second example is:

 Table[SeriesCoefficient[EllipticTheta[3, 0, x^7] EllipticTheta[3, 0, x^9] + 1/2 EllipticTheta[2, 0, x^7] EllipticTheta[2, 0, x^9], {x, 0, 2 n}], {n, 0, 10}]
 (* {1, 0, 2, 0, 0, 0, 0, 0, 4, 2, 0} *)

Under versions 7, 10.2 or 11 is all correct, but in 11.3 I get

 {1 + 2 (x^7)^(1/4) (x^9)^(1/4), 0, 0, 0, 0, 0, 0, 2 (x^7)^(1/4) (x^9)^(1/4), 4, 2 (x^7)^(1/4) (x^9)^(1/4), 0}

Bug ?

Answer 1

[Too long for a comment]

The question is based on a faulty premise, to the effect that the earlier Series-related results were correct. I'll show otherwise here.

Start with the function and subtract off the first few terms of the series exmapnsion. These are the ones that are "well behaved" in the sense that they are legitimate polynomial terms.

ee = (1/2) x^(-1/
       4) (EllipticTheta[2, 0, x^1] EllipticTheta[3, 0, x^21] - 
       EllipticTheta[3, 0, x^1] EllipticTheta[2, 0, 
         x^21])/(QPochhammer[x^3] QPochhammer[x^7]) - (1 + x^2 + x^3);

Take a series expansion to order 6. Notice we now have fractional powers. This is in version 11.3.

ss = Normal[Series[ee, {x, 0, 6}]]

(* Out[245]= x^5 + 3 x^6 - (x^21)^(1/4)/x^(1/4) - 
 2 x^(3/4) (x^21)^(1/4) - x^(11/4) (x^21)^(1/4) - 
 4 x^(15/4) (x^21)^(1/4) - 2 x^(23/4) (x^21)^(1/4) *)

In earlier versions the result for ss was simply x^6.

Now check the extent to which series approximation agrees numerically with the actual function, at eight points equally spaced on a small circle around the origin. We divide actual value by series approximation and expect to see a list of values all near 1.

ee/ss /. x -> .1*Exp[2*I*Pi*Range[0, 7]/8]

(* Out[251]= {1.12691142042, 0.995108313946 + 0.0050674357517 I, 
 0.99562715648 + 0.00461368642745 I, 
 0.932949663626 + 0.0548646912519 I, 
 1.00566321794 - 7.39707727244*10^-12 I, 
 0.932949663626 - 0.0548646912519 I, 
 0.99562715648 - 0.00461368642745 I, 
 0.995108313946 - 0.0050674357517 I} *)

Earlier version result is not so nice.

Out[4]= {1.11111, 3.07986 - 16.0451 I, 12.9841 - 7.91498 I, 
                                                 -10
    0.929903 + 0.0614126 I, -15.0789 + 1.10911 10    I, 
     0.929903 - 0.0614126 I, 12.9841 + 7.91498 I, 3.07986 + 16.0451 I}