Skip to content

Commit 432d1ff

Browse files
steppisteppi
committed
BUG: special: Fix behavior of gamma and gammasgn at poles of the Gamma function (#21827)
Fixes the behavior of `special.gamma` and `special.gammasgn` at poles of the `gamma` function. As described in #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 #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.
1 parent 16a1772 commit 432d1ff

File tree

10 files changed

+109
-81
lines changed

10 files changed

+109
-81
lines changed

‎scipy/special/xsf/cephes/beta.h

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,7 @@
5353
#include "../config.h"
5454
#include "const.h"
5555
#include "gamma.h"
56+
#include "rgamma.h"
5657

5758
namespace xsf {
5859
namespace cephes {
@@ -155,17 +156,18 @@ namespace cephes {
155156
return (sign * std::exp(y));
156157
}
157158

158-
y = Gamma(y);
159+
y = rgamma(y);
159160
a = Gamma(a);
160161
b = Gamma(b);
161-
if (y == 0.0)
162+
if (std::isinf(y)) {
162163
goto overflow;
164+
}
163165

164-
if (std::abs(std::abs(a) - std::abs(y)) > std::abs(std::abs(b) - std::abs(y))) {
165-
y = b / y;
166+
if (std::abs(std::abs(a*y) - 1.0) > std::abs(std::abs(b*y) - 1.0)) {
167+
y = b * y;
166168
y *= a;
167169
} else {
168-
y = a / y;
170+
y = a * y;
169171
y *= b;
170172
}
171173

@@ -228,20 +230,20 @@ namespace cephes {
228230
return (y);
229231
}
230232

231-
y = Gamma(y);
233+
y = rgamma(y);
232234
a = Gamma(a);
233235
b = Gamma(b);
234-
if (y == 0.0) {
236+
if (std::isinf(y)) {
235237
over:
236238
set_error("lbeta", SF_ERROR_OVERFLOW, NULL);
237239
return (sign * std::numeric_limits<double>::infinity());
238240
}
239241

240-
if (std::abs(std::abs(a) - std::abs(y)) > std::abs(std::abs(b) - std::abs(y))) {
241-
y = b / y;
242+
if (std::abs(std::abs(a*y) - 1.0) > std::abs(std::abs(b*y) - 1.0)) {
243+
y = b * y;
242244
y *= a;
243245
} else {
244-
y = a / y;
246+
y = a * y;
245247
y *= b;
246248
}
247249

‎scipy/special/xsf/cephes/expn.h

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,7 @@
6565
#include "../config.h"
6666
#include "../error.h"
6767
#include "const.h"
68-
#include "gamma.h"
68+
#include "rgamma.h"
6969
#include "polevl.h"
7070

7171
namespace xsf {
@@ -252,7 +252,7 @@ namespace cephes {
252252
k = xk;
253253
t = n;
254254
r = n - 1;
255-
ans = (std::pow(z, r) * psi / Gamma(t)) - ans;
255+
ans = (std::pow(z, r) * psi * rgamma(t)) - ans;
256256
return ans;
257257
}
258258

‎scipy/special/xsf/cephes/gamma.h

Lines changed: 37 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -149,17 +149,30 @@ namespace cephes {
149149
int sgngam = 1;
150150

151151
if (!std::isfinite(x)) {
152-
return x;
152+
if (x > 0) {
153+
// gamma(+inf) = +inf
154+
return x;
155+
}
156+
// gamma(NaN) and gamma(-inf) both should equal NaN.
157+
return std::numeric_limits<double>::quiet_NaN();
153158
}
159+
160+
if (x == 0) {
161+
/* For pole at zero, value depends on sign of zero.
162+
* +inf when approaching from right, -inf when approaching
163+
* from left. */
164+
return std::copysign(std::numeric_limits<double>::infinity(), x);
165+
}
166+
154167
q = std::abs(x);
155168

156169
if (q > 33.0) {
157170
if (x < 0.0) {
158171
p = std::floor(q);
159172
if (p == q) {
160-
gamnan:
161-
set_error("Gamma", SF_ERROR_OVERFLOW, NULL);
162-
return (std::numeric_limits<double>::infinity());
173+
// x is a negative integer. This is a pole.
174+
set_error("Gamma", SF_ERROR_SINGULAR, NULL);
175+
return (std::numeric_limits<double>::quiet_NaN());
163176
}
164177
i = p;
165178
if ((i & 1) == 0) {
@@ -215,7 +228,9 @@ namespace cephes {
215228

216229
small:
217230
if (x == 0.0) {
218-
goto gamnan;
231+
/* For this to have happened, x must have started as a negative integer. */
232+
set_error("Gamma", SF_ERROR_SINGULAR, NULL);
233+
return (std::numeric_limits<double>::quiet_NaN());
219234
} else
220235
return (z / ((1.0 + 0.5772156649015329 * x) * x));
221236
}
@@ -360,16 +375,23 @@ namespace cephes {
360375
}
361376
if (x > 0) {
362377
return 1.0;
363-
} else {
364-
fx = std::floor(x);
365-
if (x - fx == 0.0) {
366-
return 0.0;
367-
} else if (static_cast<int>(fx) % 2) {
368-
return -1.0;
369-
} else {
370-
return 1.0;
371-
}
372-
}
378+
}
379+
if (x == 0) {
380+
return std::copysign(1.0, x);
381+
}
382+
if (std::isinf(x)) {
383+
// x > 0 case handled, so x must be negative infinity.
384+
return std::numeric_limits<double>::quiet_NaN();
385+
}
386+
fx = std::floor(x);
387+
if (x - fx == 0.0) {
388+
return std::numeric_limits<double>::quiet_NaN();
389+
}
390+
// sign of gamma for x in (-n, -n+1) for positive integer n is (-1)^n.
391+
if (static_cast<int>(fx) % 2) {
392+
return -1.0;
393+
}
394+
return 1.0;
373395
}
374396

375397
} // namespace cephes

‎scipy/special/xsf/cephes/hyp2f1.h

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,7 @@
7575

7676
#include "const.h"
7777
#include "gamma.h"
78+
#include "rgamma.h"
7879
#include "psi.h"
7980

8081
namespace xsf {
@@ -252,9 +253,9 @@ namespace cephes {
252253
/* sum for t = 0 */
253254
y = xsf::cephes::psi(1.0) + xsf::cephes::psi(1.0 + e) - xsf::cephes::psi(a + d1) -
254255
xsf::cephes::psi(b + d1) - ax;
255-
y /= xsf::cephes::Gamma(e + 1.0);
256+
y *= xsf::cephes::rgamma(e + 1.0);
256257

257-
p = (a + d1) * (b + d1) * s / xsf::cephes::Gamma(e + 2.0); /* Poch for t=1 */
258+
p = (a + d1) * (b + d1) * s * xsf::cephes::rgamma(e + 2.0); /* Poch for t=1 */
258259
t = 1.0;
259260
do {
260261
r = xsf::cephes::psi(1.0 + t) + xsf::cephes::psi(1.0 + t + e) -
@@ -292,10 +293,10 @@ namespace cephes {
292293
}
293294
nosum:
294295
p = xsf::cephes::Gamma(c);
295-
y1 *= xsf::cephes::Gamma(e) * p /
296-
(xsf::cephes::Gamma(a + d1) * xsf::cephes::Gamma(b + d1));
296+
y1 *= xsf::cephes::Gamma(e) * p *
297+
(xsf::cephes::rgamma(a + d1) * xsf::cephes::rgamma(b + d1));
297298

298-
y *= p / (xsf::cephes::Gamma(a + d2) * xsf::cephes::Gamma(b + d2));
299+
y *= p * (xsf::cephes::rgamma(a + d2) * xsf::cephes::rgamma(b + d2));
299300
if ((aid & 1) != 0)
300301
y = -y;
301302

@@ -493,8 +494,8 @@ namespace cephes {
493494
p *= std::pow(-x, -a);
494495
q *= std::pow(-x, -b);
495496
t1 = Gamma(c);
496-
s = t1 * Gamma(b - a) / (Gamma(b) * Gamma(c - a));
497-
y = t1 * Gamma(a - b) / (Gamma(a) * Gamma(c - b));
497+
s = t1 * Gamma(b - a) * (rgamma(b) * rgamma(c - a));
498+
y = t1 * Gamma(a - b) * (rgamma(a) * rgamma(c - b));
498499
return s * p + y * q;
499500
} else if (x < -1.0) {
500501
if (std::abs(a) < std::abs(b)) {
@@ -532,7 +533,7 @@ namespace cephes {
532533
}
533534
if (d <= 0.0)
534535
goto hypdiv;
535-
y = Gamma(c) * Gamma(d) / (Gamma(p) * Gamma(r));
536+
y = Gamma(c) * Gamma(d) * (rgamma(p) * rgamma(r));
536537
goto hypdon;
537538
}
538539
if (d <= -1.0)

‎scipy/special/xsf/cephes/hyperg.h

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,7 @@
6767

6868
#include "const.h"
6969
#include "gamma.h"
70+
#include "rgamma.h"
7071

7172
namespace xsf {
7273
namespace cephes {
@@ -203,14 +204,14 @@ namespace cephes {
203204

204205
h1 = hyp2f0(a, a - b + 1, -1.0 / x, 1, &err1);
205206

206-
temp = std::exp(u) / xsf::cephes::Gamma(b - a);
207+
temp = std::exp(u) * xsf::cephes::rgamma(b - a);
207208
h1 *= temp;
208209
err1 *= temp;
209210

210211
h2 = hyp2f0(b - a, 1.0 - a, 1.0 / x, 2, &err2);
211212

212213
if (a < 0)
213-
temp = std::exp(t) / xsf::cephes::Gamma(a);
214+
temp = std::exp(t) * xsf::cephes::rgamma(a);
214215
else
215216
temp = std::exp(t - xsf::cephes::lgam(a));
216217

‎scipy/special/xsf/cephes/incbet.h

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,21 @@ namespace cephes {
7474

7575
namespace detail {
7676

77+
/* Compute (u * v) * w, disabling optimizations for gcc on 32 bit systems.
78+
* Used below in incbet_pseries to prevent aggressive optimizations from
79+
* degrading accuracy.
80+
*/
81+
#if defined(__GNUC__) && defined(__i386__)
82+
#pragma GCC push_options
83+
#pragma GCC optimize("00")
84+
#endif
85+
XSF_HOST_DEVICE inline double triple_product(double u, double v, double w) {
86+
return (u * v) * w;
87+
}
88+
#if defined(__GNUC__) && defined(__i386__)
89+
#pragma GCC pop_options
90+
#endif
91+
7792
constexpr double incbet_big = 4.503599627370496e15;
7893
constexpr double incbet_biginv = 2.22044604925031308085e-16;
7994

@@ -104,7 +119,7 @@ namespace cephes {
104119
u = a * std::log(x);
105120
if ((a + b) < MAXGAM && std::abs(u) < MAXLOG) {
106121
t = 1.0 / beta(a, b);
107-
s = s * t * std::pow(x, a);
122+
s = triple_product(s, t, std::pow(x, a)); // (s * t) * std::pow(x, a)
108123
} else {
109124
t = -lbeta(a, b) + u + std::log(s);
110125
if (t < MINLOG) {

‎scipy/special/xsf/cephes/jv.h

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,7 @@
5555

5656
#include "airy.h"
5757
#include "cbrt.h"
58-
#include "gamma.h"
58+
#include "rgamma.h"
5959
#include "j0.h"
6060
#include "j1.h"
6161
#include "polevl.h"
@@ -241,7 +241,7 @@ namespace cephes {
241241
t = std::frexp(0.5 * x, &ex);
242242
ex = ex * n;
243243
if ((ex > -1023) && (ex < 1023) && (n > 0.0) && (n < (MAXGAM - 1.0))) {
244-
t = std::pow(0.5 * x, n) / xsf::cephes::Gamma(n + 1.0);
244+
t = std::pow(0.5 * x, n) * xsf::cephes::rgamma(n + 1.0);
245245
y *= t;
246246
} else {
247247
t = n * std::log(0.5 * x) - lgam_sgn(n + 1.0, &sgngam);
@@ -592,13 +592,13 @@ namespace cephes {
592592

593593
if (x == 0 && n < 0 && !nint) {
594594
set_error("Jv", SF_ERROR_OVERFLOW, NULL);
595-
return std::numeric_limits<double>::infinity() / Gamma(n + 1);
595+
return std::numeric_limits<double>::infinity() * rgamma(n + 1);
596596
}
597597

598598
y = std::abs(x);
599599

600600
if (y * y < std::abs(n + 1) * detail::MACHEP) {
601-
return std::pow(0.5 * x, n) / Gamma(n + 1);
601+
return std::pow(0.5 * x, n) * rgamma(n + 1);
602602
}
603603

604604
k = 3.6 * std::sqrt(y);

‎scipy/special/xsf/cephes/rgamma.h

Lines changed: 13 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -73,35 +73,21 @@ namespace cephes {
7373

7474
XSF_HOST_DEVICE double rgamma(double x) {
7575
double w, y, z;
76-
int sign;
7776

78-
if (x > 34.84425627277176174) {
79-
return std::exp(-lgam(x));
80-
}
81-
if (x < -34.034) {
82-
w = -x;
83-
z = sinpi(w);
84-
if (z == 0.0) {
85-
return 0.0;
86-
}
87-
if (z < 0.0) {
88-
sign = 1;
89-
z = -z;
90-
} else {
91-
sign = -1;
92-
}
77+
if (x == 0) {
78+
// This case is separate from below to get correct sign for zero.
79+
return x;
80+
}
81+
82+
if (x < 0 && x == std::floor(x)) {
83+
// Gamma poles.
84+
return 0.0;
85+
}
86+
87+
if (std::abs(x) > 4.0) {
88+
return 1.0 / Gamma(x);
89+
}
9390

94-
y = std::log(w * z) - std::log(M_PI) + lgam(w);
95-
if (y < -detail::MAXLOG) {
96-
set_error("rgamma", SF_ERROR_UNDERFLOW, NULL);
97-
return (sign * 0.0);
98-
}
99-
if (y > detail::MAXLOG) {
100-
set_error("rgamma", SF_ERROR_OVERFLOW, NULL);
101-
return (sign * std::numeric_limits<double>::infinity());
102-
}
103-
return (sign * std::exp(y));
104-
}
10591
z = 1.0;
10692
w = x;
10793

‎scipy/special/xsf/cephes/struve.h

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@
8888

8989
#include "dd_real.h"
9090
#include "gamma.h"
91+
#include "rgamma.h"
9192
#include "scipy_iv.h"
9293

9394
namespace xsf {
@@ -306,7 +307,7 @@ namespace cephes {
306307
if (v < -1) {
307308
return xsf::cephes::gammasgn(v + 1.5) * std::numeric_limits<double>::infinity();
308309
} else if (v == -1) {
309-
return 2 / std::sqrt(M_PI) / xsf::cephes::Gamma(0.5);
310+
return 2 / std::sqrt(M_PI) * xsf::cephes::rgamma(0.5);
310311
} else {
311312
return 0;
312313
}

0 commit comments

Comments
 (0)