Computes the natural logarithm of #gamma(double) . This method is more numerically stable than taking the log of the result of Γ(z). Special Values: Double#POSITIVE_INFINITY returned if z ≤ 0 Double#NaN returned if z = {@link Double#NEGATIVE_INFINITY
(double z)
| 134 | * @return Log(Γ(z)) |
| 135 | */ |
| 136 | public static double lnGamma(double z) |
| 137 | { |
| 138 | if(z == Double.NEGATIVE_INFINITY) |
| 139 | return Double.NaN; |
| 140 | else if(z == Double.POSITIVE_INFINITY) |
| 141 | return z; |
| 142 | else if(z <= 0) |
| 143 | return Double.POSITIVE_INFINITY; |
| 144 | |
| 145 | /* |
| 146 | * Lanczos approximation for the log of the gamma function, with |error| < 10^-15 for all z > 0 (Almost full double precision) |
| 147 | */ |
| 148 | |
| 149 | int j; |
| 150 | double x, tmp, y, ser = 0.999999999999997092; |
| 151 | double[] c = new double[] |
| 152 | { |
| 153 | 57.1562356658629235, -59.5979603554754912, |
| 154 | 14.1360979747417471, -0.491913816097620199, .339946499848118887e-4, |
| 155 | .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3, |
| 156 | -.210264441724104883e-3, .217439618115212643e-3, -.164318106536763890e-3, |
| 157 | .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5 |
| 158 | }; |
| 159 | |
| 160 | y = x = z; |
| 161 | tmp = x+671.0/128.0; |
| 162 | tmp = (x+0.5)*log(tmp)-tmp; |
| 163 | for (j = 0; j < 14; j++) |
| 164 | { |
| 165 | y++; |
| 166 | ser += c[j] / y; |
| 167 | } |
| 168 | |
| 169 | return tmp+log(2.5066282746310005*ser/x); |
| 170 | } |
| 171 | |
| 172 | /** |
| 173 | * Positive zero of the digamma function |
no test coverage detected