diff options
author | nobu <nobu@b2dd03c8-39d4-4d8f-98ff-823fe69b080e> | 2017-05-12 16:01:17 +0000 |
---|---|---|
committer | nobu <nobu@b2dd03c8-39d4-4d8f-98ff-823fe69b080e> | 2017-05-12 16:01:17 +0000 |
commit | 05ea264611b5b8f73d75c5676478d401420419ce (patch) | |
tree | 43f17cacf910c4ecc51cc175902f3c5f5fd09607 /missing/tgamma.c | |
parent | 2c0baa97a99e5fd89693e740fa5ef3ed7f23e6d7 (diff) | |
download | ruby-05ea264611b5b8f73d75c5676478d401420419ce.tar.gz |
tgamma.c: unify versions with/without lgamma_r
git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@58689 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
Diffstat (limited to 'missing/tgamma.c')
-rw-r--r-- | missing/tgamma.c | 36 |
1 files changed, 10 insertions, 26 deletions
diff --git a/missing/tgamma.c b/missing/tgamma.c index 5e306fbb43..9324e19ea1 100644 --- a/missing/tgamma.c +++ b/missing/tgamma.c @@ -13,30 +13,7 @@ reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten #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,15 +45,16 @@ 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 (x < 0) { - int sign; static double zero = 0.0; double i, f; f = modf(-x, &i); @@ -84,9 +62,15 @@ double tgamma(double x) /* Gamma function */ 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 +} |