[388] | 1 | /*
|
---|
| 2 | "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
|
---|
| 3 | "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
|
---|
| 4 | "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
|
---|
| 5 |
|
---|
| 6 | approximation method:
|
---|
| 7 |
|
---|
| 8 | (x - 0.5) S(x)
|
---|
| 9 | Gamma(x) = (x + g - 0.5) * ----------------
|
---|
| 10 | exp(x + g - 0.5)
|
---|
| 11 |
|
---|
| 12 | with
|
---|
| 13 | a1 a2 a3 aN
|
---|
| 14 | S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
|
---|
| 15 | x + 1 x + 2 x + 3 x + N
|
---|
| 16 |
|
---|
| 17 | with a0, a1, a2, a3,.. aN constants which depend on g.
|
---|
| 18 |
|
---|
| 19 | for x < 0 the following reflection formula is used:
|
---|
| 20 |
|
---|
| 21 | Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
|
---|
| 22 |
|
---|
| 23 | most ideas and constants are from boost and python
|
---|
| 24 | */
|
---|
| 25 | #include "libm.h"
|
---|
| 26 |
|
---|
| 27 | static const double pi = 3.141592653589793238462643383279502884;
|
---|
| 28 |
|
---|
| 29 | /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
|
---|
| 30 | static double sinpi(double x)
|
---|
| 31 | {
|
---|
| 32 | int n;
|
---|
| 33 |
|
---|
| 34 | /* argument reduction: x = |x| mod 2 */
|
---|
| 35 | /* spurious inexact when x is odd int */
|
---|
| 36 | x = x * 0.5;
|
---|
| 37 | x = 2 * (x - floor(x));
|
---|
| 38 |
|
---|
| 39 | /* reduce x into [-.25,.25] */
|
---|
| 40 | n = 4 * x;
|
---|
| 41 | n = (n+1)/2;
|
---|
| 42 | x -= n * 0.5;
|
---|
| 43 |
|
---|
| 44 | x *= pi;
|
---|
| 45 | switch (n) {
|
---|
| 46 | default: /* case 4 */
|
---|
| 47 | case 0:
|
---|
| 48 | return __sin(x, 0, 0);
|
---|
| 49 | case 1:
|
---|
| 50 | return __cos(x, 0);
|
---|
| 51 | case 2:
|
---|
| 52 | return __sin(-x, 0, 0);
|
---|
| 53 | case 3:
|
---|
| 54 | return -__cos(x, 0);
|
---|
| 55 | }
|
---|
| 56 | }
|
---|
| 57 |
|
---|
| 58 | #define N 12
|
---|
| 59 | //static const double g = 6.024680040776729583740234375;
|
---|
| 60 | static const double gmhalf = 5.524680040776729583740234375;
|
---|
| 61 | static const double Snum[N+1] = {
|
---|
| 62 | 23531376880.410759688572007674451636754734846804940,
|
---|
| 63 | 42919803642.649098768957899047001988850926355848959,
|
---|
| 64 | 35711959237.355668049440185451547166705960488635843,
|
---|
| 65 | 17921034426.037209699919755754458931112671403265390,
|
---|
| 66 | 6039542586.3520280050642916443072979210699388420708,
|
---|
| 67 | 1439720407.3117216736632230727949123939715485786772,
|
---|
| 68 | 248874557.86205415651146038641322942321632125127801,
|
---|
| 69 | 31426415.585400194380614231628318205362874684987640,
|
---|
| 70 | 2876370.6289353724412254090516208496135991145378768,
|
---|
| 71 | 186056.26539522349504029498971604569928220784236328,
|
---|
| 72 | 8071.6720023658162106380029022722506138218516325024,
|
---|
| 73 | 210.82427775157934587250973392071336271166969580291,
|
---|
| 74 | 2.5066282746310002701649081771338373386264310793408,
|
---|
| 75 | };
|
---|
| 76 | static const double Sden[N+1] = {
|
---|
| 77 | 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535,
|
---|
| 78 | 2637558, 357423, 32670, 1925, 66, 1,
|
---|
| 79 | };
|
---|
| 80 | /* n! for small integer n */
|
---|
| 81 | static const double fact[] = {
|
---|
| 82 | 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0,
|
---|
| 83 | 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0,
|
---|
| 84 | 355687428096000.0, 6402373705728000.0, 121645100408832000.0,
|
---|
| 85 | 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
|
---|
| 86 | };
|
---|
| 87 |
|
---|
| 88 | /* S(x) rational function for positive x */
|
---|
| 89 | static double S(double x)
|
---|
| 90 | {
|
---|
| 91 | double_t num = 0, den = 0;
|
---|
| 92 | int i;
|
---|
| 93 |
|
---|
| 94 | /* to avoid overflow handle large x differently */
|
---|
| 95 | if (x < 8)
|
---|
| 96 | for (i = N; i >= 0; i--) {
|
---|
| 97 | num = num * x + Snum[i];
|
---|
| 98 | den = den * x + Sden[i];
|
---|
| 99 | }
|
---|
| 100 | else
|
---|
| 101 | for (i = 0; i <= N; i++) {
|
---|
| 102 | num = num / x + Snum[i];
|
---|
| 103 | den = den / x + Sden[i];
|
---|
| 104 | }
|
---|
| 105 | return num/den;
|
---|
| 106 | }
|
---|
| 107 |
|
---|
| 108 | double tgamma(double x)
|
---|
| 109 | {
|
---|
| 110 | union {double f; uint64_t i;} u = {x};
|
---|
| 111 | double absx, y;
|
---|
| 112 | double_t dy, z, r;
|
---|
| 113 | uint32_t ix = u.i>>32 & 0x7fffffff;
|
---|
| 114 | int sign = u.i>>63;
|
---|
| 115 |
|
---|
| 116 | /* special cases */
|
---|
| 117 | if (ix >= 0x7ff00000)
|
---|
| 118 | /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
|
---|
| 119 | return x + INFINITY;
|
---|
| 120 | if (ix < (0x3ff-54)<<20)
|
---|
| 121 | /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
|
---|
| 122 | return 1/x;
|
---|
| 123 |
|
---|
| 124 | /* integer arguments */
|
---|
| 125 | /* raise inexact when non-integer */
|
---|
| 126 | if (x == floor(x)) {
|
---|
| 127 | if (sign)
|
---|
| 128 | return 0/0.0;
|
---|
| 129 | if (x <= sizeof fact/sizeof *fact)
|
---|
| 130 | return fact[(int)x - 1];
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | /* x >= 172: tgamma(x)=inf with overflow */
|
---|
| 134 | /* x =< -184: tgamma(x)=+-0 with underflow */
|
---|
| 135 | if (ix >= 0x40670000) { /* |x| >= 184 */
|
---|
| 136 | if (sign) {
|
---|
| 137 | FORCE_EVAL((float)(0x1p-126/x));
|
---|
| 138 | if (floor(x) * 0.5 == floor(x * 0.5))
|
---|
| 139 | return 0;
|
---|
| 140 | return -0.0;
|
---|
| 141 | }
|
---|
| 142 | x *= 0x1p1023;
|
---|
| 143 | return x;
|
---|
| 144 | }
|
---|
| 145 |
|
---|
| 146 | absx = sign ? -x : x;
|
---|
| 147 |
|
---|
| 148 | /* handle the error of x + g - 0.5 */
|
---|
| 149 | y = absx + gmhalf;
|
---|
| 150 | if (absx > gmhalf) {
|
---|
| 151 | dy = y - absx;
|
---|
| 152 | dy -= gmhalf;
|
---|
| 153 | } else {
|
---|
| 154 | dy = y - gmhalf;
|
---|
| 155 | dy -= absx;
|
---|
| 156 | }
|
---|
| 157 |
|
---|
| 158 | z = absx - 0.5;
|
---|
| 159 | r = S(absx) * exp(-y);
|
---|
| 160 | if (x < 0) {
|
---|
| 161 | /* reflection formula for negative x */
|
---|
| 162 | /* sinpi(absx) is not 0, integers are already handled */
|
---|
| 163 | r = -pi / (sinpi(absx) * absx * r);
|
---|
| 164 | dy = -dy;
|
---|
| 165 | z = -z;
|
---|
| 166 | }
|
---|
| 167 | r += dy * (gmhalf+0.5) * r / y;
|
---|
| 168 | z = pow(y, 0.5*z);
|
---|
| 169 | y = r * z * z;
|
---|
| 170 | return y;
|
---|
| 171 | }
|
---|
| 172 |
|
---|
| 173 | #if 0
|
---|
| 174 | double __lgamma_r(double x, int *sign)
|
---|
| 175 | {
|
---|
| 176 | double r, absx;
|
---|
| 177 |
|
---|
| 178 | *sign = 1;
|
---|
| 179 |
|
---|
| 180 | /* special cases */
|
---|
| 181 | if (!isfinite(x))
|
---|
| 182 | /* lgamma(nan)=nan, lgamma(+-inf)=inf */
|
---|
| 183 | return x*x;
|
---|
| 184 |
|
---|
| 185 | /* integer arguments */
|
---|
| 186 | if (x == floor(x) && x <= 2) {
|
---|
| 187 | /* n <= 0: lgamma(n)=inf with divbyzero */
|
---|
| 188 | /* n == 1,2: lgamma(n)=0 */
|
---|
| 189 | if (x <= 0)
|
---|
| 190 | return 1/0.0;
|
---|
| 191 | return 0;
|
---|
| 192 | }
|
---|
| 193 |
|
---|
| 194 | absx = fabs(x);
|
---|
| 195 |
|
---|
| 196 | /* lgamma(x) ~ -log(|x|) for tiny |x| */
|
---|
| 197 | if (absx < 0x1p-54) {
|
---|
| 198 | *sign = 1 - 2*!!signbit(x);
|
---|
| 199 | return -log(absx);
|
---|
| 200 | }
|
---|
| 201 |
|
---|
| 202 | /* use tgamma for smaller |x| */
|
---|
| 203 | if (absx < 128) {
|
---|
| 204 | x = tgamma(x);
|
---|
| 205 | *sign = 1 - 2*!!signbit(x);
|
---|
| 206 | return log(fabs(x));
|
---|
| 207 | }
|
---|
| 208 |
|
---|
| 209 | /* second term (log(S)-g) could be more precise here.. */
|
---|
| 210 | /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
|
---|
| 211 | r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5));
|
---|
| 212 | if (x < 0) {
|
---|
| 213 | /* reflection formula for negative x */
|
---|
| 214 | x = sinpi(absx);
|
---|
| 215 | *sign = 2*!!signbit(x) - 1;
|
---|
| 216 | r = log(pi/(fabs(x)*absx)) - r;
|
---|
| 217 | }
|
---|
| 218 | return r;
|
---|
| 219 | }
|
---|
| 220 |
|
---|
| 221 | weak_alias(__lgamma_r, lgamma_r);
|
---|
| 222 | #endif
|
---|