| OLD | NEW |
| 1 /* | 1 /* |
| 2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964) | 2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964) |
| 3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001) | 3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001) |
| 4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004) | 4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004) |
| 5 | 5 |
| 6 approximation method: | 6 approximation method: |
| 7 | 7 |
| 8 (x - 0.5) S(x) | 8 (x - 0.5) S(x) |
| 9 Gamma(x) = (x + g - 0.5) * ---------------- | 9 Gamma(x) = (x + g - 0.5) * ---------------- |
| 10 exp(x + g - 0.5) | 10 exp(x + g - 0.5) |
| (...skipping 172 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
| 183 dy = -dy; | 183 dy = -dy; |
| 184 z = -z; | 184 z = -z; |
| 185 } | 185 } |
| 186 r += dy * (gmhalf + 0.5) * r / y; | 186 r += dy * (gmhalf + 0.5) * r / y; |
| 187 z = pow(y, 0.5 * z); | 187 z = pow(y, 0.5 * z); |
| 188 y = r * z * z; | 188 y = r * z * z; |
| 189 return y; | 189 return y; |
| 190 } | 190 } |
| 191 | 191 |
| 192 #if 0 | 192 #if 0 |
| 193 double __lgamma_r(double x, int *sign) | 193 double __lgamma_r(double x, int* sign) { |
| 194 { | 194 double r, absx; |
| 195 » double r, absx; | |
| 196 | 195 |
| 197 » *sign = 1; | 196 *sign = 1; |
| 198 | 197 |
| 199 » /* special cases */ | 198 /* special cases */ |
| 200 » if (!isfinite(x)) | 199 if (!isfinite(x)) |
| 201 » » /* lgamma(nan)=nan, lgamma(+-inf)=inf */ | 200 /* lgamma(nan)=nan, lgamma(+-inf)=inf */ |
| 202 » » return x*x; | 201 return x * x; |
| 203 | 202 |
| 204 » /* integer arguments */ | 203 /* integer arguments */ |
| 205 » if (x == floor(x) && x <= 2) { | 204 if (x == floor(x) && x <= 2) { |
| 206 » » /* n <= 0: lgamma(n)=inf with divbyzero */ | 205 /* n <= 0: lgamma(n)=inf with divbyzero */ |
| 207 » » /* n == 1,2: lgamma(n)=0 */ | 206 /* n == 1,2: lgamma(n)=0 */ |
| 208 » » if (x <= 0) | 207 if (x <= 0) |
| 209 » » » return 1/0.0; | 208 return 1 / 0.0; |
| 210 » » return 0; | 209 return 0; |
| 211 » } | 210 } |
| 212 | 211 |
| 213 » absx = fabs(x); | 212 absx = fabs(x); |
| 214 | 213 |
| 215 » /* lgamma(x) ~ -log(|x|) for tiny |x| */ | 214 /* lgamma(x) ~ -log(|x|) for tiny |x| */ |
| 216 » if (absx < 0x1p-54) { | 215 if (absx < 0x1p-54) { |
| 217 » » *sign = 1 - 2*!!signbit(x); | 216 *sign = 1 - 2 * !!signbit(x); |
| 218 » » return -log(absx); | 217 return -log(absx); |
| 219 » } | 218 } |
| 220 | 219 |
| 221 » /* use tgamma for smaller |x| */ | 220 /* use tgamma for smaller |x| */ |
| 222 » if (absx < 128) { | 221 if (absx < 128) { |
| 223 » » x = tgamma(x); | 222 x = tgamma(x); |
| 224 » » *sign = 1 - 2*!!signbit(x); | 223 *sign = 1 - 2 * !!signbit(x); |
| 225 » » return log(fabs(x)); | 224 return log(fabs(x)); |
| 226 » } | 225 } |
| 227 | 226 |
| 228 » /* second term (log(S)-g) could be more precise here.. */ | 227 /* second term (log(S)-g) could be more precise here.. */ |
| 229 » /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */ | 228 /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */ |
| 230 » r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5)); | 229 r = (absx - 0.5) * (log(absx + gmhalf) - 1) + (log(S(absx)) - (gmhalf + 0.5)); |
| 231 » if (x < 0) { | 230 if (x < 0) { |
| 232 » » /* reflection formula for negative x */ | 231 /* reflection formula for negative x */ |
| 233 » » x = sinpi(absx); | 232 x = sinpi(absx); |
| 234 » » *sign = 2*!!signbit(x) - 1; | 233 *sign = 2 * !!signbit(x) - 1; |
| 235 » » r = log(pi/(fabs(x)*absx)) - r; | 234 r = log(pi / (fabs(x) * absx)) - r; |
| 236 » } | 235 } |
| 237 » return r; | 236 return r; |
| 238 } | 237 } |
| 239 | 238 |
| 240 weak_alias(__lgamma_r, lgamma_r); | 239 weak_alias(__lgamma_r, lgamma_r); |
| 241 #endif | 240 #endif |
| OLD | NEW |