summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authornobu <nobu@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2017-05-12 16:01:17 +0000
committernobu <nobu@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2017-05-12 16:01:17 +0000
commit05ea264611b5b8f73d75c5676478d401420419ce (patch)
tree43f17cacf910c4ecc51cc175902f3c5f5fd09607
parent2c0baa97a99e5fd89693e740fa5ef3ed7f23e6d7 (diff)
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
-rw-r--r--missing/tgamma.c36
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
+}