Repository: mpmath/mpmath (python)

Issue:
Unexpected Result from laguerre with inf input
mpmath version: 1.1.0

I'm no expert on this function, but it appears that the expected behavior for `laguerre(1, a, b)` is to return `1 + a - b`. SciPy agrees with this. However, when `b` is `inf`, the return appears to consistently be `1 + a` when one might expect a negative infinity. This is what SciPy and GSL do in this situation. Interestingly, the fp version of laguerre takes a different route entirely, throwing a `NoConvergence` exception.

```
Python 3.6.8 (default, Oct  7 2019, 12:59:55) 
[GCC 8.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from mpmath import *
>>> from scipy import special
>>> mp.laguerre(1, 2, 5)
mpf('-2.0')
>>> mp.laguerre(1, 2, mp.inf)
mpf('3.0')
>>> mp.laguerre(1, 7.2, mp.inf)
mpf('8.1999999999999993')
>>> fp.laguerre(1, 7.2, fp.inf)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.6/dist-packages/mpmath/ctx_fp.py", line 57, in f_wrapped
    return f(ctx, *args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/mpmath/functions/orthogonal.py", line 348, in laguerre
    return ctx.hypercomb(h, [a], **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/mpmath/functions/hypergeometric.py", line 127, in hypercomb
    [ctx.rgamma(b) for b in beta_s] + \
  File "/usr/local/lib/python3.6/dist-packages/mpmath/functions/hypergeometric.py", line 222, in hyper
    if   q == 1: return ctx._hyp1f1(a_s, b_s, z, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/mpmath/functions/hypergeometric.py", line 342, in _hyp1f1
    v = ctx.hypsum(1, 1, (atype, btype), [a, b], z, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/mpmath/ctx_fp.py", line 213, in hypsum
    raise ctx.NoConvergence
mpmath.libmp.libhyper.NoConvergence
>>> special.eval_genlaguerre(1, 2, fp.inf)
-inf
>>> special.eval_genlaguerre(1, 7.2, fp.inf)
-inf
```

Pull request:
Special case in ctx.hypsum for infinite z
Closes #511

Patch:
diff --git a/mpmath/ctx_fp.py b/mpmath/ctx_fp.py
index 79cb6ba..9e931cf 100644
--- a/mpmath/ctx_fp.py
+++ b/mpmath/ctx_fp.py
@@ -196,16 +196,40 @@ def _is_real_type(ctx, z):
     def _is_complex_type(ctx, z):
         return isinstance(z, complex)
 
-    def hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs):
-        coeffs = list(coeffs)
+    def hypsum(ctx, p, q, flags, coeffs, z, maxterms=6000, **kwargs):
+        for i, c in enumerate(coeffs[p:], start=p):
+            if flags[i] == 'Z':
+                if c <= 0:
+                    ok = False
+                    for ii, cc in enumerate(coeffs[:p]):
+                        # Note: c <= cc or c < cc, depending on convention
+                        if flags[ii] == 'Z' and cc <= 0 and c <= cc:
+                            ok = True
+                    if not ok:
+                        raise ZeroDivisionError("pole in hypergeometric series")
         num = range(p)
         den = range(p,p+q)
+        if ctx.isinf(z):
+            n = max(((n, c) for n, c in enumerate(coeffs[:p])
+                     if flags[n] == 'Z' and c < 0), default=(-1, 0),
+                    key=lambda x: x[1])[0]
+            if n >= 0:
+                n = -coeffs[n]
+                t = z**n
+                for k in range(n):
+                    for i in num: t *= (coeffs[i]+k)
+                    for i in den: t /= (coeffs[i]+k)
+                    t /= (k+1)
+                return t
         tol = ctx.eps
         s = t = 1.0
         k = 0
         while 1:
             for i in num: t *= (coeffs[i]+k)
-            for i in den: t /= (coeffs[i]+k)
+            try:
+                for i in den: t /= (coeffs[i]+k)
+            except ZeroDivisionError:
+                raise NotImplementedError
             k += 1; t /= k; t *= z; s += t
             if abs(t) < tol:
                 return s
diff --git a/mpmath/ctx_mp.py b/mpmath/ctx_mp.py
index d815dde..d2fde03 100644
--- a/mpmath/ctx_mp.py
+++ b/mpmath/ctx_mp.py
@@ -645,6 +645,30 @@ def hypsum(ctx, p, q, flags, coeffs, z, accurate_small=True, **kwargs):
         elif hasattr(z, "_mpc_"):
             key = p, q, flags, 'C'
             v = z._mpc_
+        for i, c in enumerate(coeffs[p:], start=p):
+            if flags[i] == 'Z':
+                if c <= 0:
+                    ok = False
+                    for ii, cc in enumerate(coeffs[:p]):
+                        # Note: c <= cc or c < cc, depending on convention
+                        if flags[ii] == 'Z' and cc <= 0 and c <= cc:
+                            ok = True
+                    if not ok:
+                        raise ZeroDivisionError("pole in hypergeometric series")
+        num = range(p)
+        den = range(p,p+q)
+        if ctx.isinf(z):
+            n = max(((n, c) for n, c in enumerate(coeffs[:p])
+                     if flags[n] == 'Z' and c < 0), default=(-1, 0),
+                    key=lambda x: x[1])[0]
+            if n >= 0:
+                n = -coeffs[n]
+                t = z**n
+                for k in range(n):
+                    for i in num: t *= (coeffs[i]+k)
+                    for i in den: t /= (coeffs[i]+k)
+                    t /= (k+1)
+                return t
         if key not in ctx.hyp_summators:
             ctx.hyp_summators[key] = libmp.make_hyp_summator(key)[1]
         summator = ctx.hyp_summators[key]
@@ -659,14 +683,6 @@ def hypsum(ctx, p, q, flags, coeffs, z, accurate_small=True, **kwargs):
         max_total_jump = 0
         for i, c in enumerate(coeffs):
             if flags[i] == 'Z':
-                if i >= p and c <= 0:
-                    ok = False
-                    for ii, cc in enumerate(coeffs[:p]):
-                        # Note: c <= cc or c < cc, depending on convention
-                        if flags[ii] == 'Z' and cc <= 0 and c <= cc:
-                            ok = True
-                    if not ok:
-                        raise ZeroDivisionError("pole in hypergeometric series")
                 continue
             n, d = ctx.nint_distance(c)
             n = -int(n)


Test patch:
diff --git a/mpmath/tests/test_functions2.py b/mpmath/tests/test_functions2.py
index c57ec95..1176fb4 100644
--- a/mpmath/tests/test_functions2.py
+++ b/mpmath/tests/test_functions2.py
@@ -165,6 +165,8 @@ def test_hyper_misc():
     mp.dps = 25
     v = mpc('1.2282306665029814734863026', '-0.1225033830118305184672133')
     assert hyper([(3,4),2+j,1],[1,5,j/3],mpf(1)/5+j/8).ae(v)
+    pytest.raises(ZeroDivisionError, lambda: mp.hyper([1, 2, -2], [-1, 3], 1.1))
+    pytest.raises(ZeroDivisionError, lambda: fp.hyper([1, 2, -2], [-1, 3], 1.1))
 
 def test_elliptic_integrals():
     assert ellipk(0).ae(pi/2)
@@ -528,6 +530,7 @@ def test_hyper_2f1():
 def test_hyper_2f1_hard():
     # Singular cases
     assert hyp2f1(2,-1,-1,3).ae(7)
+    pytest.raises(NotImplementedError, lambda: fp.hyp2f1(2,-1,-1,3))
     assert hyp2f1(2,-1,-1,3,eliminate_all=True).ae(0.25)
     assert hyp2f1(2,-2,-2,3).ae(34)
     assert hyp2f1(2,-2,-2,3,eliminate_all=True).ae(0.25)
@@ -2412,3 +2415,8 @@ def test_issue_505():
 
 def test_issue_653():
     pytest.raises(ValueError, lambda: zeta(2, -2))
+
+def test_issue_511():
+    assert mp.laguerre(1, 2, mp.inf) == -mp.inf
+    assert mp.laguerre(1, 7.2, mp.inf) == -mp.inf
+    assert fp.laguerre(1, 7.2, fp.inf) == -fp.inf