Avoid spurious overflow in fp gammaprod (#1150)

* Avoid spurious overflow in fp gammaprod

fp.binomial(1100, 1) raised OverflowError even though the result is
exactly 1100: gammaprod multiplies the individual gamma values, and
gamma(1101) exceeds the double range although the quotient does not.

Fall back to evaluating the regular numerator/denominator in log space
when a term overflows, so the quotient is computed whenever it is
representable. Genuinely out-of-range results, such as
fp.binomial(1100, 550), still raise OverflowError.

Closes #493

mpmath/functions/factorials.py

@@ -31,8 +31,19 @@ def gammaprod(ctx, a, b, _infsign=False):
             i = poles_num.pop()
             j = poles_den.pop()
             p *= (-1)**(i+j) * ctx.gamma(1-j) / ctx.gamma(1-i)
-        for x in regular_num: p *= ctx.gamma(x)
-        for x in regular_den: p /= ctx.gamma(x)
+        try:
+            q = ctx.one
+            for x in regular_num: q *= ctx.gamma(x)
+            for x in regular_den: q /= ctx.gamma(x)
+        except OverflowError:
+            # In the fp context an individual gamma value can exceed the
+            # double range even when the quotient is representable, e.g.
+            # binomial(1100, 1). Evaluate the regular part in log space.
+            s = ctx.zero
+            for x in regular_num: s += ctx.loggamma(x)
+            for x in regular_den: s -= ctx.loggamma(x)
+            q = ctx.exp(s)
+        p *= q
     finally:
         ctx.prec = orig
     return +p

mpmath/tests/test_fp.py

@@ -1811,3 +1811,12 @@ def test_issue_491():
 def test_issue_521():
     assert fp.ff(1, -fp.inf) == 0.0
     assert fp.isnan(fp.ff(1, fp.inf))
+
+def test_issue_493():
+    assert ae(fp.binomial(1100, 1), 1100.0)
+    assert ae(fp.binomial(1100, 1099), 1100.0)
+    assert fp.binomial(1100, 0) == 1.0
+    assert ae(fp.rf(1100, 1), 1100.0)
+    assert ae(fp.beta(1100, 1), 1/1100)
+    assert ae(fp.binomial(5, 2), 10.0)
+    pytest.raises(OverflowError, lambda: fp.binomial(1100, 550))