diff options
| -rw-r--r-- | ChangeLog | 5 | ||||
| -rw-r--r-- | missing/lgamma_r.c | 9 |
2 files changed, 10 insertions, 4 deletions
@@ -1,3 +1,8 @@ +Sat Feb 9 11:09:26 2008 Tanaka Akira <akr@fsij.org> + + * missing/xlgamma_r.c (lgamma_r): return HUGE_VAL for non-positive + integers. + Sat Feb 9 10:03:07 2008 Tanaka Akira <akr@fsij.org> * string.c (rb_str_new4): copy encoding from orig, instead of shared diff --git a/missing/lgamma_r.c b/missing/lgamma_r.c index d4402c67b..ab6c28294 100644 --- a/missing/lgamma_r.c +++ b/missing/lgamma_r.c @@ -12,6 +12,7 @@ reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten gamma.c -- Gamma function ***********************************************************/ #include <math.h> +#include <errno.h> #define PI 3.14159265358979324 /* $\pi$ */ #define LOG_2PI 1.83787706640934548 /* $\log 2\pi$ */ #define LOG_PI 1.14472988584940017 /* $\log_e \pi$ */ @@ -47,13 +48,13 @@ loggamma(double x) /* the natural logarithm of the Gamma function. */ double lgamma_r(double x, int *signp) { - if (x < 0) { + if (x <= 0) { double i, f, s; f = modf(-x, &i); - if (f == 0.0) { - static const double zero = 0.0; + if (f == 0.0) { /* pole error */ *signp = 1; - return 1.0/zero; + errno = ERANGE; + return HUGE_VAL; } *signp = (fmod(i, 2.0) != 0.0) ? 1 : -1; s = sin(PI * f); |