Computes the value of the digamma function, Ψ(x), which is the derivative of #lnGamma(double) . This method may return: Double#NaN for zero and negative integer values (would be complex infinity) This method should be accurate to an absolu
(double x)
| 190 | * @return the value of Ψ(x) |
| 191 | */ |
| 192 | public static double digamma(double x) |
| 193 | { |
| 194 | if(x == 0) |
| 195 | return Double.NaN;//complex infinity |
| 196 | else if(x < 0)//digamma(1-x) == digamma(x)+pi/tan(pi*x), to make x positive |
| 197 | { |
| 198 | if(Math.rint(x) == x) |
| 199 | return Double.NaN;//the zeros are complex infinity |
| 200 | return digamma(1-x)-PI/tan(PI*x); |
| 201 | } |
| 202 | else if(x < 2)//shift the value into [2, Inf] |
| 203 | return digamma(x+1) - 1/x; |
| 204 | double approx= log(x);//rel error of 10^-11 for x >= 100 10^-13 for x >= 250, near machine precision at x >= 500 |
| 205 | double approxS = -1/(2*x); |
| 206 | double approxSS = -1/(12*x*x); |
| 207 | |
| 208 | if(x <= 7) |
| 209 | return (x-digammaPosZero)*hornerPolyR(digamma_p_2_7, x)/hornerPolyR(digamma_q_2_7, x); |
| 210 | else if(x <= 70) |
| 211 | return approx + (approxS + (approxSS + hornerPolyR(digamma_p_7_70, x)/hornerPolyR(digamma_q_7_70, x))); |
| 212 | if(x < 500) |
| 213 | return approx + (approxS + (approxSS + hornerPolyR(digamma_adj_p, x)/hornerPolyR(digamma_adj_q, x))); |
| 214 | else |
| 215 | return approx; |
| 216 | } |
| 217 | |
| 218 | /** |
| 219 | * The upper polynomial for adjustment to the digamma function approximation |