* Compute the natural log of Gamma(x), accurate to 10 decimal places. * * This implementation is based on: * * Pike, M.C., I.D. Hill (1966) Algorithm 291: Logarithm of Gamma function * [S14]. Communications of the ACM 9(9):684. */
| 7 | * [S14]. Communications of the ACM 9(9):684. |
| 8 | */ |
| 9 | static inline double |
| 10 | ln_gamma(double x) { |
| 11 | double f, z; |
| 12 | |
| 13 | assert(x > 0.0); |
| 14 | |
| 15 | if (x < 7.0) { |
| 16 | f = 1.0; |
| 17 | z = x; |
| 18 | while (z < 7.0) { |
| 19 | f *= z; |
| 20 | z += 1.0; |
| 21 | } |
| 22 | x = z; |
| 23 | f = -log(f); |
| 24 | } else { |
| 25 | f = 0.0; |
| 26 | } |
| 27 | |
| 28 | z = 1.0 / (x * x); |
| 29 | |
| 30 | return f + (x-0.5) * log(x) - x + 0.918938533204673 + |
| 31 | (((-0.000595238095238 * z + 0.000793650793651) * z - |
| 32 | 0.002777777777778) * z + 0.083333333333333) / x; |
| 33 | } |
| 34 | |
| 35 | /* |
| 36 | * Compute the incomplete Gamma ratio for [0..x], where p is the shape |
no outgoing calls
no test coverage detected