Incorrect limit involving elliptic functions

[oscarbenjamin]

This is from Sage:
https://groups.google.com/g/sage-support/c/dDe6JLMOgbM/m/mQLMehDFAAAJ?utm_medium=email&utm_source=footer

e = ((1 - 3*x**2)*elliptic_e(4*x/(x**2 + 2*x + 1))**2/2 - (x**2 - 2*x + 1)*elliptic_e(4*x/(x**2 + 2*x + 1))*elliptic_k(4*x/(x**2 + 2*x + 1))/2)/(pi**2*x**8 - 2*pi**2*x**7 - pi**2*x**6 + 4*pi**2*x**5 - pi**2*x**4 - 2*pi**2*x**3 + pi**2*x**2)

That gives:

In [9]: e.subs(x, 0.01).n()
Out[9]: -0.125070333715528

In [10]: e.limit(x, 0)
Out[10]: -5/8

The answer should be -1/8.

[anutosh491]

The issues here come out of gruntz (which tells me there is something wrong with the series expansion of the numerator)

>>> gruntz(e, x, 0)
-5/8
>>> num
(1 - 3*x**2)*elliptic_e(4*x/(x**2 + 2*x + 1))**2/2 - (x**2 - 2*x + 1)*elliptic_e(4*x/(x**2 + 2*x + 1))*elliptic_k(4*x/(x**2 + 2*x + 1))/2
>>> num.series(x, 0, 3)
-5*pi**2*x**2/8 + O(x**3)

I think this is where the -5/8 originates from. If someone is debugging this, try the following

1. I think we don't have as_leading_term for elliptic_e & elliptic_k implemented as of now. So try getting this in place first
2. Then look into what's wrong with the series expansion. (it faulty term might either be generated by elliptic_e or elliptic_k)

[haru-44]

Here is a smaller example where a similar phenomenon occurs.

>>> x = symbols('x')
>>> t = 4*x/(x+1)
>>> e = ((x+1)*elliptic_e(t) + (x-1)*elliptic_k(t))/x**2
>>> print(e.limit(x, 0).evalf())
-4.71238898038469
>>> print(e.subs(x, 0.0001).evalf())
-1.57111057497161

We use nseries to evaluate the limit.

>>> print(e.nseries().subs(x, 0).evalf())
-4.71238898038469

So how about the following example?

>>> elliptic_e(t).nseries()

$$\frac{\pi}{2} - \frac{\pi x}{2} - \frac{3 \pi x^{2}}{8} - \frac{5 \pi x^{3}}{8} - \frac{175 \pi x^{4}}{128} - \frac{441 \pi x^{5}}{128} + O\left(x^{6}\right)$$

This differs from the wolfram results.

$$\frac{\pi}{2} - \frac{\pi x}{2} + \frac{\pi x^{2}}{8} - \frac{3 \pi x^{3}}{8} - \frac{15 \pi x^{4}}{128} + O\left(x^{5}\right)$$

[oscarbenjamin]

It's a bug in hyper.nseries:

In [105]: f = hyper((-S(1)/2, S(1)/2), (1,), 4*x/(x + 1))

In [106]: f.nseries()
Out[106]: 
                  2            3             4             5           
      x        3⋅x          5⋅x         175⋅x         441⋅x        ⎛ 6⎞
1 - ───── - ────────── - ────────── - ─────────── - ─────────── + O⎝x ⎠
    x + 1            2            3             4             5        
            4⋅(x + 1)    4⋅(x + 1)    64⋅(x + 1)    64⋅(x + 1) 

The method does not seem to do anything to handle the argument being anything other than a symbol.

This change gives correct answers for all cases listed above:

diff --git a/sympy/functions/special/hyper.py b/sympy/functions/special/hyper.py
index 277a1dda48..3ddb3cef19 100644
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -282,7 +282,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):
         for i in range(n):
             num = Mul(*[RisingFactorial(a, i) for a in ap])
             den = Mul(*[RisingFactorial(b, i) for b in bq])
-            terms.append(((num/den) * (arg**i)) / factorial(i))
+            terms.append(((num/den) * (arg**i).nseries()) / factorial(i))
 
         return (Add(*terms) + Order(x**n,x))

I'm not sure what the proper way to write this code is.

This also works:

diff --git a/sympy/functions/special/hyper.py b/sympy/functions/special/hyper.py
index 277a1dda48..7b75611249 100644
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -265,27 +265,6 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             return S.One
         return super()._eval_as_leading_term(x, logx=logx, cdir=cdir)
 
-    def _eval_nseries(self, x, n, logx, cdir=0):
-
-        from sympy.series.order import Order
-
-        arg = self.args[2]
-        x0 = arg.limit(x, 0)
-        ap = self.args[0]
-        bq = self.args[1]
-
-        if x0 != 0:
-            return super()._eval_nseries(x, n, logx)
-
-        terms = []
-
-        for i in range(n):
-            num = Mul(*[RisingFactorial(a, i) for a in ap])
-            den = Mul(*[RisingFactorial(b, i) for b in bq])
-            terms.append(((num/den) * (arg**i)) / factorial(i))
-
-        return (Add(*terms) + Order(x**n,x))
-
     @property
     def argument(self):
         """ Argument of the hypergeometric function. """

Maybe _eval_nseries should not be defined by hyper.

[haru-44]

hyper._eval_nseries implements expression

$$\,{}_pF_q(a_1,\ldots,a_p;b_1,\ldots,b_q;x) = \sum_{n=0}^\infty \frac{(a_1)_n\cdots(a_p)_n}{(b_1)_n\cdots(b_q)_n} \, \frac {x^n} {n!}$$

from Wikipedia.

However, it does not assume a case where the argument is not $x$ but $x/(x+1)$. So it seems to me that I should change (arg**i) to (arg**i).nseries(), but I can't be sure that this would work for all cases.

On the other hand, Function._eval_nseries computes series using diff even if _eval_nseries is not implemented. Then it seems that we don't have to implement _eval_nseries, but in fact there is a difference in speed.

%%timeit
hyper((-S(1)/2, S(1)/2), (1,), x)._eval_nseries(x, n=50, logx=None, cdir=0)
880 µs ± 40.6 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%%timeit
hyperexpand(Function._eval_nseries(hyper((-S(1)/2, S(1)/2), (1,), x), x, n=50, logx=None, cdir=0))
6.41 ms ± 260 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

The safest way would be to evaluate with hyper._eval_nseries only when the argument is $x$, and with Function._eval_nseries otherwise.

[oscarbenjamin]

So it seems to me that I should change (arg**i) to (arg**i).nseries(), but I can't be sure that this would work for all cases.

It would probably be best to compute arg.nseries() and check that it is an ordinary power series with no singularities, negative powers, log terms etc. Then the formula could be implemented using the terms from arg.nseries().

In any case the bug that gives incorrect results should be fixed so if there is any doubt about what to do it would be better to begin by deleting the method that gives incorrect results.

[haru-44]

Yes, it is. I think we can delete it.