diff options
Diffstat (limited to 'missing/tgamma.c')
| -rw-r--r-- | missing/tgamma.c | 44 |
1 files changed, 17 insertions, 27 deletions
diff --git a/missing/tgamma.c b/missing/tgamma.c index 5e306fbb43..82d614d755 100644 --- a/missing/tgamma.c +++ b/missing/tgamma.c @@ -9,34 +9,12 @@ reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten /*********************************************************** gamma.c -- Gamma function ***********************************************************/ -#include "ruby/config.h" +#include "ruby/internal/config.h" +#include "ruby/missing.h" #include <math.h> #include <errno.h> -#ifdef HAVE_LGAMMA_R - -double tgamma(double x) -{ - int sign; - double d; - if (x == 0.0) { /* Pole Error */ - errno = ERANGE; - return 1/x < 0 ? -HUGE_VAL : HUGE_VAL; - } - if (x < 0) { - static double zero = 0.0; - double i, f; - f = modf(-x, &i); - if (f == 0.0) { /* Domain Error */ - errno = EDOM; - return zero/zero; - } - } - d = lgamma_r(x, &sign); - return sign * exp(d); -} - -#else +#ifndef HAVE_LGAMMA_R #include <errno.h> #define PI 3.14159265358979324 /* $\pi$ */ @@ -68,25 +46,37 @@ loggamma(double x) /* the natural logarithm of the Gamma function. */ + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x + 0.5 * LOG_2PI - log(v) - x + (x - 0.5) * log(x); } +#endif double tgamma(double x) /* Gamma function */ { + int sign; if (x == 0.0) { /* Pole Error */ errno = ERANGE; return 1/x < 0 ? -HUGE_VAL : HUGE_VAL; } + if (isinf(x)) { + if (x < 0) goto domain_error; + return x; + } if (x < 0) { - int sign; static double zero = 0.0; double i, f; f = modf(-x, &i); if (f == 0.0) { /* Domain Error */ + domain_error: errno = EDOM; return zero/zero; } +#ifndef HAVE_LGAMMA_R sign = (fmod(i, 2.0) != 0.0) ? 1 : -1; return sign * PI / (sin(PI * f) * exp(loggamma(1 - x))); +#endif } +#ifndef HAVE_LGAMMA_R return exp(loggamma(x)); -} +#else + x = lgamma_r(x, &sign); + return sign * exp(x); #endif +} |