# Is it possible to circumvent a bug inside SeriesCoefficient?

Bug persisting through 13.1.

As far as I can tell, there seems to be a bug in SeriesCoefficient:

In[4]:=
f[z_] := (-z + z Sqrt[1 - 4 z^2])/(2 (-1 + z));
SeriesCoefficient[f[z], {z, 0, 10}]
SeriesCoefficient[f[z], {z, 0, n}] /. n -> 10

Out[5]= 9

Out[6]= 93/256


Note that both outputs should be the same, but apparently the application general in n makes a mistake (9 is the correct result). In fact, I know that the coefficients grow exponentially in n but we get:

In[9]:=
SeriesCoefficient[f[z], {z, 0, n},
Assumptions :> Mod[n, 2] == 0 && n >= 2] // FullSimplify

Out[9]= 1/2 (1 - Gamma[n/2]/(Sqrt[\[Pi]] Gamma[(1 + n)/2]))


This quantity,

$$\qquad\displaystyle \frac{1}{2} \cdot \biggl(1 - \frac{\Gamma\bigl(\frac{n}{2}\bigr)}{\sqrt{\pi} \cdot \Gamma\bigr(\frac{n+1}{2}\bigr)} \biggr)$$,

converges to 1/2 from below (odd n support the same series), so it's not even close to being correct.

Obviously, I can't tell what exactly goes wrong inside SeriesCoefficient.

Since I have an interest in the general form of the coefficient (of the expansions of this and other functions), is there a way to circumvent the bug inside (or the application of SeriesCoefficient) so that it works more reliably? Less power would be acceptable in trade.

• Mathematica 11.0 outputs $\frac{1}{2} \cdot (1 - i \sqrt{3})$ for all $n \geq 1$ -- still wrong, but different. Apr 2, 2017 at 16:26

As

l = CoefficientList[Series[f[z], {z, 0, 20}], z]
(*
{0, 0, 0, 1, 1, 2, 2, 4, 4, 9, 9, 23, 23, 65, 65, 197, 197, 626, 626, 2056, 2056}
*)


we can calculate only the even terms.

You could do for example:

k = FindSequenceFunction[l[[2 ;; ;; 2]]]

Chop@N@k[5]  (* The tenth coefficient *)
(*
9
*)


Moreover:

FullSimplify[k[n], Element[k, Integers] && k > 0]
(*
(2-2 I Sqrt@3-(4^n Gamma[-1/2+n] Hypergeometric2F1Regularized[1,-1/2+n,1+n, 4])/Sqrt@π)/4
*)


$$k(n) = \frac{1}{4} \left(-\frac{4^n \, _2\tilde{F}_1\left(1,n-\frac{1}{2};n+1;4\right) \Gamma \left(n-\frac{1}{2}\right)}{\sqrt{\pi }}-2 i \sqrt{3}+2\right)$$

• Nice, never knew of these functions. One concern: FindSequenceFunction tries to find a function for which the finite list of coefficients matches, right? Isn't that inherently unsave? May 9, 2014 at 9:19
• @Raphael Of course! But you may try to "confirm" the result by testing on higher coefficients. I think it's pretty safe in this case (but no, I wont bet my money on it :) ) May 9, 2014 at 9:34
• I guess more can't be accepted at this point. (It is possible to prove the result by induction, but that's a mess I'd like to avoid.) Since I kind of have to bet my money on the result, I guess I'll have to wait for the bugfix. Thanks! May 9, 2014 at 16:59

The basic idea is to replace the 1/(1-z) in f[z] with its series and multiply with the series of the complementary factor of f[z]. The series coefficient of f[z] is then a sum of the series coefficients of the complementary factor. It turns out Mathematica can handle this. (Er, I think...see end.) Like @belisarius, we find a formula for every other coefficient, since the SeriesCoefficient[] below is zero whenever n-k is even.

sc0 = Sum[
SeriesCoefficient[
(z - z Sqrt[1 - 4 z^2])/2*z^k,
{z, 0, n},
Assumptions -> k ∈ Integers && n >= k >= 0] /. {k ->
2 j, n -> 2 m + 1},
{j, 0, m}] // FullSimplify[#, m ∈ Integers && m > 0] &
(*
1/2 (1 - (-4)^m *
Binomial[1/2, m] Hypergeometric2F1[1, -m, 3/2 - m, 1/4])
*)


Test sc0 against SeriesCoefficient of f[z]:

Table[sc0 /. m -> Max[Ceiling[n/2 - 1], 0], {n, 0, 20}]
test1 = Table[ (* the original f[z] at explicit degrees *)
SeriesCoefficient[
(-z + z Sqrt[1 - 4 z^2])/(2 (-1 + z)),
{z, 0, n}],
{n, 0, 20}];
test0 == test1
(*
{0, 0, 0, 1, 1, 2, 2, 4, 4, 9, 9, 23, 23,
65, 65, 197, 197, 626, 626, 2056, 2056}
True
*)


Updated explanation. It's an induction proof that sc0 is equivalent to @belisarius's answer bel (the index is adjusted to match sc0). Checking the induction step is apparently susceptible to Fudd's First Law of Opposition: Push hard enough on the simplification and the difference of the differences will eventual collapse.

bel = (2 - 2 I Sqrt@3 -
(4^n Gamma[-1/2 + n] *
Hypergeometric2F1Regularized[1, -1/2 + n, 1 + n, 4])/Sqrt@π)/4 /.
n -> m + 1;
sc0 - bel /. m -> 0 // Simplify
sc0 - bel /. {{m -> m + 1}, {m -> m}} // Apply@Subtract //
FunctionExpand // FullSimplify[#, m ∈ Integers && m > 0] & //
FunctionExpand // FullSimplify[#, m ∈ Integers && m > 0] &
(*
0  <-- base step
0  <-- induction step
*)


Er, I think....The first trial had the first coefficient wrong (my mistake). A tiny little mod that should have fixed it resulted in all complex-number coefficients, not integers. Going back to what I thought I started with resulted in different complex coefficients. Having turned back down the wrong path, I turned around again and decided to forge ahead. After trying a very clever :) trick that revealed several related bugs in Sum, I managed to construct a double sum that Mathematica could do accurately. Since some time had passed since this morning, I though let's just try the original idea again. And violà! you have the answer above.