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 223