BUG: special: Fix behavior of gamma and gammasgn at poles of the Gamma function (scipy#21827)

Fixes the behavior of `special.gamma` and `special.gammasgn` at poles of the `gamma`
function. As described in scipy#21810, we previously had `special.gamma(n) == np.inf` for
every non-positive integer `n`, and bizarrely `special.gamma(-np.inf) == -np.inf`.
Similarly, as described in scipy#19710, we've had `gammasgn(n) == 0.0` for every
non-positive integer `n` and even more bizarrely, `gammasgn(-np.inf) == 1.0`.

The sign of $\Gamma(x)$ at each pole depends on the direction in which $x$ approaches
the pole, so it doesn't make sense to return $+\infty$ at all poles. Accordingly, the
C99 standard recommends the following behavior for the gamma function in special cases.

1. If the argument is $\pm 0$, $\pm \infty$ is returned.
2. If the argument is a negative integer, NaN is returned.
3. If the argument is $-\infty$, NaN is returned.
4. If the argument is $+\infty$, $+\infty$ is returned
5. If the argument is NaN, NaN is returned.

This updates `special.gamma` to follow this standard. Note that signed zeros offer the
ability to specify the sign of infinity of `gamma(x)` at `x == 0`. This also updates
`gammasgn` to be consistent with the above.

Additionally, improve the accuracy of the `rgamma` implementation, and use that as a
replacement where we relied on the previous behaviour; document this in the docstring.

Author: steppi

scipy/special/_special_ufuncs_docs.cpp

@@ -1572,10 +1572,35 @@ const char *gamma_doc = R"(
     which, combined with the fact that :math:`\Gamma(1) = 1`, implies
     the above identity for :math:`z = n`.
 
+    The gamma function has poles at non-negative integers and the sign
+    of infinity as z approaches each pole depends upon the direction in
+    which the pole is approached. For this reason, the consistent thing
+    is for gamma(z) to return NaN at negative integers, and to return
+    -inf when x = -0.0 and +inf when x = 0.0, using the signbit of zero
+    to signify the direction in which the origin is being approached. This
+    is for instance what is recommended for the gamma function in annex F
+    entry 9.5.4 of the Iso C 99 standard [isoc99]_.
+
+    Prior to SciPy version 1.15, ``scipy.special.gamma(z)`` returned ``+inf``
+    at each pole. This was fixed in version 1.15, but with the following
+    consequence. Expressions where gamma appears in the denominator
+    such as
+
+    ``gamma(u) * gamma(v) / (gamma(w) * gamma(x))``
+
+    no longer evaluate to 0 if the numerator is well defined but there is a
+    pole in the denominator. Instead such expressions evaluate to NaN. We
+    recommend instead using the function `rgamma` for the reciprocal gamma
+    function in such cases. The above expression could for instance be written
+    as
+
+    ``gamma(u) * gamma(v) * (rgamma(w) * rgamma(x))``
+
     References
     ----------
     .. [dlmf] NIST Digital Library of Mathematical Functions
               https://dlmf.nist.gov/5.2#E1
+    .. [isoc99] https://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf
 
     Examples
     --------

scipy/special/tests/test_basic.py

@@ -398,8 +398,13 @@ def test_gammaln(self):
         cephes.gammaln(10)
 
     def test_gammasgn(self):
-        vals = np.array([-4, -3.5, -2.3, 1, 4.2], np.float64)
-        assert_array_equal(cephes.gammasgn(vals), np.sign(cephes.rgamma(vals)))
+        vals = np.array(
+            [-np.inf, -4, -3.5, -2.3, -0.0, 0.0, 1, 4.2, np.inf], np.float64
+        )
+        reference = np.array(
+            [np.nan, np.nan, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0], np.float64
+        )
+        assert_array_equal(cephes.gammasgn(vals), reference)
 
     def test_gdtr(self):
         assert_equal(cephes.gdtr(1,1,0),0.0)
@@ -2715,9 +2720,26 @@ def test_rgamma(self):
         assert_almost_equal(rgam,rlgam,8)
 
     def test_infinity(self):
-        assert_(np.isinf(special.gamma(-1)))
         assert_equal(special.rgamma(-1), 0)
 
+    @pytest.mark.parametrize(
+        "x,expected",
+        [
+            # infinities
+            ([-np.inf, np.inf], [np.nan, np.inf]),
+            # negative and positive zero
+            ([-0.0, 0.0], [-np.inf, np.inf]),
+            # small poles
+            (range(-32, 0), np.full(32, np.nan)),
+            # medium sized poles
+            (range(-1024, -32, 99), np.full(11, np.nan)),
+            # large pole
+            ([-4.141512231792294e+16], [np.nan]),
+        ]
+    )
+    def test_poles(self, x, expected):
+        assert_array_equal(special.gamma(x), expected)
+
 
 class TestHankel:
 
@@ -3690,7 +3712,7 @@ def test_pbdv_seq(self):
     def test_pbdv_points(self):
         # simple case
         eta = np.linspace(-10, 10, 5)
-        z = 2**(eta/2)*np.sqrt(np.pi)/special.gamma(.5-.5*eta)
+        z = 2**(eta/2)*np.sqrt(np.pi)*special.rgamma(.5-.5*eta)
         assert_allclose(special.pbdv(eta, 0.)[0], z, rtol=1e-14, atol=1e-14)
 
         # some points

scipy/special/tests/test_hyp2f1.py

@@ -2488,6 +2488,27 @@ def test_region6(self, hyp2f1_test_case):
         )
         assert_allclose(hyp2f1(a, b, c, z), expected, rtol=rtol)
 
+
+    @pytest.mark.parametrize(
+        "hyp2f1_test_case",
+        [
+            # Broke when fixing gamma pole behavior in gh-21827
+            pytest.param(
+                Hyp2f1TestCase(
+                    a=1.3,
+                    b=-0.2,
+                    c=0.3,
+                    z=-2.1,
+                    expected=1.8202169687521206,
+                    rtol=5e-15,
+                ),
+            ),
+        ]
+    )
+    def test_miscellaneous(self, hyp2f1_test_case ):
+        a, b, c, z, expected, rtol = hyp2f1_test_case
+        assert_allclose(hyp2f1(a, b, c, z), expected, rtol=rtol)
+
     @pytest.mark.slow
     @check_version(mpmath, "1.0.0")
     def test_test_hyp2f1(self):

scipy/special/xsf/cephes/beta.h

@@ -53,6 +53,7 @@
 #include "../config.h"
 #include "const.h"
 #include "gamma.h"
+#include "rgamma.h"
 
 namespace xsf {
 namespace cephes {
@@ -155,17 +156,18 @@ namespace cephes {
             return (sign * std::exp(y));
         }
 
-        y = Gamma(y);
+        y = rgamma(y);
         a = Gamma(a);
         b = Gamma(b);
-        if (y == 0.0)
+        if (std::isinf(y)) {
             goto overflow;
+	}
 
-        if (std::abs(std::abs(a) - std::abs(y)) > std::abs(std::abs(b) - std::abs(y))) {
-            y = b / y;
+        if (std::abs(std::abs(a*y) - 1.0) > std::abs(std::abs(b*y) - 1.0)) {
+            y = b * y;
             y *= a;
         } else {
-            y = a / y;
+            y = a * y;
             y *= b;
         }
 
@@ -228,20 +230,20 @@ namespace cephes {
             return (y);
         }
 
-        y = Gamma(y);
+        y = rgamma(y);
         a = Gamma(a);
         b = Gamma(b);
-        if (y == 0.0) {
+        if (std::isinf(y)) {
         over:
             set_error("lbeta", SF_ERROR_OVERFLOW, NULL);
             return (sign * std::numeric_limits<double>::infinity());
         }
 
-        if (std::abs(std::abs(a) - std::abs(y)) > std::abs(std::abs(b) - std::abs(y))) {
-            y = b / y;
+        if (std::abs(std::abs(a*y) - 1.0) > std::abs(std::abs(b*y) - 1.0)) {
+            y = b * y;
             y *= a;
         } else {
-            y = a / y;
+            y = a * y;
             y *= b;
         }
 

scipy/special/xsf/cephes/expn.h

@@ -65,7 +65,7 @@
 #include "../config.h"
 #include "../error.h"
 #include "const.h"
-#include "gamma.h"
+#include "rgamma.h"
 #include "polevl.h"
 
 namespace xsf {
@@ -252,7 +252,7 @@ namespace cephes {
         k = xk;
         t = n;
         r = n - 1;
-        ans = (std::pow(z, r) * psi / Gamma(t)) - ans;
+        ans = (std::pow(z, r) * psi * rgamma(t)) - ans;
         return ans;
     }
 

scipy/special/xsf/cephes/gamma.h

@@ -149,17 +149,30 @@ namespace cephes {
         int sgngam = 1;
 
         if (!std::isfinite(x)) {
-            return x;
+	    if (x > 0) {
+		// gamma(+inf) = +inf
+		return x;
+	    }
+	    // gamma(NaN) and gamma(-inf) both should equal NaN.
+            return std::numeric_limits<double>::quiet_NaN();
         }
+
+	if (x == 0) {
+	    /* For pole at zero, value depends on sign of zero.
+	     * +inf when approaching from right, -inf when approaching
+	     * from left. */
+	    return std::copysign(std::numeric_limits<double>::infinity(), x);
+	}
+
         q = std::abs(x);
 
         if (q > 33.0) {
             if (x < 0.0) {
                 p = std::floor(q);
                 if (p == q) {
-                gamnan:
-                    set_error("Gamma", SF_ERROR_OVERFLOW, NULL);
-                    return (std::numeric_limits<double>::infinity());
+		    // x is a negative integer. This is a pole.
+                    set_error("Gamma", SF_ERROR_SINGULAR, NULL);
+                    return (std::numeric_limits<double>::quiet_NaN());
                 }
                 i = p;
                 if ((i & 1) == 0) {
@@ -215,7 +228,9 @@ namespace cephes {
 
     small:
         if (x == 0.0) {
-            goto gamnan;
+	    /* For this to have happened, x must have started as a negative integer. */
+	    set_error("Gamma", SF_ERROR_SINGULAR, NULL);
+	    return (std::numeric_limits<double>::quiet_NaN());
         } else
             return (z / ((1.0 + 0.5772156649015329 * x) * x));
     }
@@ -360,16 +375,23 @@ namespace cephes {
         }
         if (x > 0) {
             return 1.0;
-        } else {
-            fx = std::floor(x);
-            if (x - fx == 0.0) {
-                return 0.0;
-            } else if (static_cast<int>(fx) % 2) {
-                return -1.0;
-            } else {
-                return 1.0;
-            }
-        }
+	}
+	if (x == 0) {
+	    return std::copysign(1.0, x);
+	}
+	if (std::isinf(x)) {
+	    // x > 0 case handled, so x must be negative infinity.
+	    return std::numeric_limits<double>::quiet_NaN();
+	}
+	fx = std::floor(x);
+	if (x - fx == 0.0) {
+	    return std::numeric_limits<double>::quiet_NaN();
+	}
+	// sign of gamma for x in (-n, -n+1) for positive integer n is (-1)^n.
+	if (static_cast<int>(fx) % 2) {
+	    return -1.0;
+	}
+	return 1.0;
     }
 
 } // namespace cephes

scipy/special/xsf/cephes/hyp2f1.h

@@ -75,6 +75,7 @@
 
 #include "const.h"
 #include "gamma.h"
+#include "rgamma.h"
 #include "psi.h"
 
 namespace xsf {
@@ -252,9 +253,9 @@ namespace cephes {
                     /* sum for t = 0 */
                     y = xsf::cephes::psi(1.0) + xsf::cephes::psi(1.0 + e) - xsf::cephes::psi(a + d1) -
                         xsf::cephes::psi(b + d1) - ax;
-                    y /= xsf::cephes::Gamma(e + 1.0);
+                    y *= xsf::cephes::rgamma(e + 1.0);
 
-                    p = (a + d1) * (b + d1) * s / xsf::cephes::Gamma(e + 2.0); /* Poch for t=1 */
+                    p = (a + d1) * (b + d1) * s * xsf::cephes::rgamma(e + 2.0); /* Poch for t=1 */
                     t = 1.0;
                     do {
                         r = xsf::cephes::psi(1.0 + t) + xsf::cephes::psi(1.0 + t + e) -
@@ -292,10 +293,10 @@ namespace cephes {
                     }
                 nosum:
                     p = xsf::cephes::Gamma(c);
-                    y1 *= xsf::cephes::Gamma(e) * p /
-                          (xsf::cephes::Gamma(a + d1) * xsf::cephes::Gamma(b + d1));
+                    y1 *= xsf::cephes::Gamma(e) * p *
+                          (xsf::cephes::rgamma(a + d1) * xsf::cephes::rgamma(b + d1));
 
-                    y *= p / (xsf::cephes::Gamma(a + d2) * xsf::cephes::Gamma(b + d2));
+                    y *= p * (xsf::cephes::rgamma(a + d2) * xsf::cephes::rgamma(b + d2));
                     if ((aid & 1) != 0)
                         y = -y;
 
@@ -493,8 +494,8 @@ namespace cephes {
             p *= std::pow(-x, -a);
             q *= std::pow(-x, -b);
             t1 = Gamma(c);
-            s = t1 * Gamma(b - a) / (Gamma(b) * Gamma(c - a));
-            y = t1 * Gamma(a - b) / (Gamma(a) * Gamma(c - b));
+            s = t1 * Gamma(b - a) * (rgamma(b) * rgamma(c - a));
+            y = t1 * Gamma(a - b) * (rgamma(a) * rgamma(c - b));
             return s * p + y * q;
         } else if (x < -1.0) {
             if (std::abs(a) < std::abs(b)) {
@@ -532,7 +533,7 @@ namespace cephes {
                 }
                 if (d <= 0.0)
                     goto hypdiv;
-                y = Gamma(c) * Gamma(d) / (Gamma(p) * Gamma(r));
+                y = Gamma(c) * Gamma(d) * (rgamma(p) * rgamma(r));
                 goto hypdon;
             }
             if (d <= -1.0)

scipy/special/xsf/cephes/hyperg.h

@@ -67,6 +67,7 @@
 
 #include "const.h"
 #include "gamma.h"
+#include "rgamma.h"
 
 namespace xsf {
 namespace cephes {
@@ -203,14 +204,14 @@ namespace cephes {
 
             h1 = hyp2f0(a, a - b + 1, -1.0 / x, 1, &err1);
 
-            temp = std::exp(u) / xsf::cephes::Gamma(b - a);
+            temp = std::exp(u) * xsf::cephes::rgamma(b - a);
             h1 *= temp;
             err1 *= temp;
 
             h2 = hyp2f0(b - a, 1.0 - a, 1.0 / x, 2, &err2);
 
             if (a < 0)
-                temp = std::exp(t) / xsf::cephes::Gamma(a);
+                temp = std::exp(t) * xsf::cephes::rgamma(a);
             else
                 temp = std::exp(t - xsf::cephes::lgam(a));
 

scipy/special/xsf/cephes/incbet.h

@@ -74,6 +74,21 @@ namespace cephes {
 
     namespace detail {
 
+	/* Compute (u * v) * w, disabling optimizations for gcc on 32 bit systems.
+	 * Used below in incbet_pseries to prevent aggressive optimizations from
+	 * degrading accuracy.
+	 */
+#if defined(__GNUC__) && defined(__i386__)
+#pragma GCC push_options
+#pragma GCC optimize("00")
+#endif
+	XSF_HOST_DEVICE inline double triple_product(double u, double v, double w) {
+	    return (u * v) * w;
+	}
+#if defined(__GNUC__) && defined(__i386__)
+#pragma GCC pop_options
+#endif
+
         constexpr double incbet_big = 4.503599627370496e15;
         constexpr double incbet_biginv = 2.22044604925031308085e-16;
 
@@ -104,7 +119,7 @@ namespace cephes {
             u = a * std::log(x);
             if ((a + b) < MAXGAM && std::abs(u) < MAXLOG) {
                 t = 1.0 / beta(a, b);
-                s = s * t * std::pow(x, a);
+                s = triple_product(s, t, std::pow(x, a)); // (s * t) * std::pow(x, a)
             } else {
                 t = -lbeta(a, b) + u + std::log(s);
                 if (t < MINLOG) {