Index: fusl/src/math/tgamma.c |
diff --git a/fusl/src/math/tgamma.c b/fusl/src/math/tgamma.c |
index 28f6e0f8327a9368c0c39bc2db4b25033cb24895..feb0616ec674db8028211c684652d0b5ae9b99f6 100644 |
--- a/fusl/src/math/tgamma.c |
+++ b/fusl/src/math/tgamma.c |
@@ -27,147 +27,166 @@ most ideas and constants are from boost and python |
static const double pi = 3.141592653589793238462643383279502884; |
/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */ |
-static double sinpi(double x) |
-{ |
- int n; |
- |
- /* argument reduction: x = |x| mod 2 */ |
- /* spurious inexact when x is odd int */ |
- x = x * 0.5; |
- x = 2 * (x - floor(x)); |
- |
- /* reduce x into [-.25,.25] */ |
- n = 4 * x; |
- n = (n+1)/2; |
- x -= n * 0.5; |
- |
- x *= pi; |
- switch (n) { |
- default: /* case 4 */ |
- case 0: |
- return __sin(x, 0, 0); |
- case 1: |
- return __cos(x, 0); |
- case 2: |
- return __sin(-x, 0, 0); |
- case 3: |
- return -__cos(x, 0); |
- } |
+static double sinpi(double x) { |
+ int n; |
+ |
+ /* argument reduction: x = |x| mod 2 */ |
+ /* spurious inexact when x is odd int */ |
+ x = x * 0.5; |
+ x = 2 * (x - floor(x)); |
+ |
+ /* reduce x into [-.25,.25] */ |
+ n = 4 * x; |
+ n = (n + 1) / 2; |
+ x -= n * 0.5; |
+ |
+ x *= pi; |
+ switch (n) { |
+ default: /* case 4 */ |
+ case 0: |
+ return __sin(x, 0, 0); |
+ case 1: |
+ return __cos(x, 0); |
+ case 2: |
+ return __sin(-x, 0, 0); |
+ case 3: |
+ return -__cos(x, 0); |
+ } |
} |
#define N 12 |
-//static const double g = 6.024680040776729583740234375; |
+// static const double g = 6.024680040776729583740234375; |
static const double gmhalf = 5.524680040776729583740234375; |
-static const double Snum[N+1] = { |
- 23531376880.410759688572007674451636754734846804940, |
- 42919803642.649098768957899047001988850926355848959, |
- 35711959237.355668049440185451547166705960488635843, |
- 17921034426.037209699919755754458931112671403265390, |
- 6039542586.3520280050642916443072979210699388420708, |
- 1439720407.3117216736632230727949123939715485786772, |
- 248874557.86205415651146038641322942321632125127801, |
- 31426415.585400194380614231628318205362874684987640, |
- 2876370.6289353724412254090516208496135991145378768, |
- 186056.26539522349504029498971604569928220784236328, |
- 8071.6720023658162106380029022722506138218516325024, |
- 210.82427775157934587250973392071336271166969580291, |
- 2.5066282746310002701649081771338373386264310793408, |
+static const double Snum[N + 1] = { |
+ 23531376880.410759688572007674451636754734846804940, |
+ 42919803642.649098768957899047001988850926355848959, |
+ 35711959237.355668049440185451547166705960488635843, |
+ 17921034426.037209699919755754458931112671403265390, |
+ 6039542586.3520280050642916443072979210699388420708, |
+ 1439720407.3117216736632230727949123939715485786772, |
+ 248874557.86205415651146038641322942321632125127801, |
+ 31426415.585400194380614231628318205362874684987640, |
+ 2876370.6289353724412254090516208496135991145378768, |
+ 186056.26539522349504029498971604569928220784236328, |
+ 8071.6720023658162106380029022722506138218516325024, |
+ 210.82427775157934587250973392071336271166969580291, |
+ 2.5066282746310002701649081771338373386264310793408, |
}; |
-static const double Sden[N+1] = { |
- 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, |
- 2637558, 357423, 32670, 1925, 66, 1, |
+static const double Sden[N + 1] = { |
+ 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535, |
+ 2637558, 357423, 32670, 1925, 66, 1, |
}; |
/* n! for small integer n */ |
static const double fact[] = { |
- 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0, |
- 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0, |
- 355687428096000.0, 6402373705728000.0, 121645100408832000.0, |
- 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0, |
+ 1, |
+ 1, |
+ 2, |
+ 6, |
+ 24, |
+ 120, |
+ 720, |
+ 5040.0, |
+ 40320.0, |
+ 362880.0, |
+ 3628800.0, |
+ 39916800.0, |
+ 479001600.0, |
+ 6227020800.0, |
+ 87178291200.0, |
+ 1307674368000.0, |
+ 20922789888000.0, |
+ 355687428096000.0, |
+ 6402373705728000.0, |
+ 121645100408832000.0, |
+ 2432902008176640000.0, |
+ 51090942171709440000.0, |
+ 1124000727777607680000.0, |
}; |
/* S(x) rational function for positive x */ |
-static double S(double x) |
-{ |
- double_t num = 0, den = 0; |
- int i; |
- |
- /* to avoid overflow handle large x differently */ |
- if (x < 8) |
- for (i = N; i >= 0; i--) { |
- num = num * x + Snum[i]; |
- den = den * x + Sden[i]; |
- } |
- else |
- for (i = 0; i <= N; i++) { |
- num = num / x + Snum[i]; |
- den = den / x + Sden[i]; |
- } |
- return num/den; |
+static double S(double x) { |
+ double_t num = 0, den = 0; |
+ int i; |
+ |
+ /* to avoid overflow handle large x differently */ |
+ if (x < 8) |
+ for (i = N; i >= 0; i--) { |
+ num = num * x + Snum[i]; |
+ den = den * x + Sden[i]; |
+ } |
+ else |
+ for (i = 0; i <= N; i++) { |
+ num = num / x + Snum[i]; |
+ den = den / x + Sden[i]; |
+ } |
+ return num / den; |
} |
-double tgamma(double x) |
-{ |
- union {double f; uint64_t i;} u = {x}; |
- double absx, y; |
- double_t dy, z, r; |
- uint32_t ix = u.i>>32 & 0x7fffffff; |
- int sign = u.i>>63; |
- |
- /* special cases */ |
- if (ix >= 0x7ff00000) |
- /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ |
- return x + INFINITY; |
- if (ix < (0x3ff-54)<<20) |
- /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ |
- return 1/x; |
- |
- /* integer arguments */ |
- /* raise inexact when non-integer */ |
- if (x == floor(x)) { |
- if (sign) |
- return 0/0.0; |
- if (x <= sizeof fact/sizeof *fact) |
- return fact[(int)x - 1]; |
- } |
- |
- /* x >= 172: tgamma(x)=inf with overflow */ |
- /* x =< -184: tgamma(x)=+-0 with underflow */ |
- if (ix >= 0x40670000) { /* |x| >= 184 */ |
- if (sign) { |
- FORCE_EVAL((float)(0x1p-126/x)); |
- if (floor(x) * 0.5 == floor(x * 0.5)) |
- return 0; |
- return -0.0; |
- } |
- x *= 0x1p1023; |
- return x; |
- } |
- |
- absx = sign ? -x : x; |
- |
- /* handle the error of x + g - 0.5 */ |
- y = absx + gmhalf; |
- if (absx > gmhalf) { |
- dy = y - absx; |
- dy -= gmhalf; |
- } else { |
- dy = y - gmhalf; |
- dy -= absx; |
- } |
- |
- z = absx - 0.5; |
- r = S(absx) * exp(-y); |
- if (x < 0) { |
- /* reflection formula for negative x */ |
- /* sinpi(absx) is not 0, integers are already handled */ |
- r = -pi / (sinpi(absx) * absx * r); |
- dy = -dy; |
- z = -z; |
- } |
- r += dy * (gmhalf+0.5) * r / y; |
- z = pow(y, 0.5*z); |
- y = r * z * z; |
- return y; |
+double tgamma(double x) { |
+ union { |
+ double f; |
+ uint64_t i; |
+ } u = {x}; |
+ double absx, y; |
+ double_t dy, z, r; |
+ uint32_t ix = u.i >> 32 & 0x7fffffff; |
+ int sign = u.i >> 63; |
+ |
+ /* special cases */ |
+ if (ix >= 0x7ff00000) |
+ /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ |
+ return x + INFINITY; |
+ if (ix < (0x3ff - 54) << 20) |
+ /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ |
+ return 1 / x; |
+ |
+ /* integer arguments */ |
+ /* raise inexact when non-integer */ |
+ if (x == floor(x)) { |
+ if (sign) |
+ return 0 / 0.0; |
+ if (x <= sizeof fact / sizeof *fact) |
+ return fact[(int)x - 1]; |
+ } |
+ |
+ /* x >= 172: tgamma(x)=inf with overflow */ |
+ /* x =< -184: tgamma(x)=+-0 with underflow */ |
+ if (ix >= 0x40670000) { /* |x| >= 184 */ |
+ if (sign) { |
+ FORCE_EVAL((float)(0x1p-126 / x)); |
+ if (floor(x) * 0.5 == floor(x * 0.5)) |
+ return 0; |
+ return -0.0; |
+ } |
+ x *= 0x1p1023; |
+ return x; |
+ } |
+ |
+ absx = sign ? -x : x; |
+ |
+ /* handle the error of x + g - 0.5 */ |
+ y = absx + gmhalf; |
+ if (absx > gmhalf) { |
+ dy = y - absx; |
+ dy -= gmhalf; |
+ } else { |
+ dy = y - gmhalf; |
+ dy -= absx; |
+ } |
+ |
+ z = absx - 0.5; |
+ r = S(absx) * exp(-y); |
+ if (x < 0) { |
+ /* reflection formula for negative x */ |
+ /* sinpi(absx) is not 0, integers are already handled */ |
+ r = -pi / (sinpi(absx) * absx * r); |
+ dy = -dy; |
+ z = -z; |
+ } |
+ r += dy * (gmhalf + 0.5) * r / y; |
+ z = pow(y, 0.5 * z); |
+ y = r * z * z; |
+ return y; |
} |
#if 0 |