The gamma function is a generalization of the factorial function. This method provides the gamma function for values from -Infinity to Infinity. Special Values: Double#NaN returned if z = 0 or z = Double#NEGATIVE_INFINITY @param z any real value @return Γ(z)
(double z)
| 65 | * @return Γ(z) |
| 66 | */ |
| 67 | public static double gamma(double z) |
| 68 | { |
| 69 | if(z == 0)//It is actualy infinity*I, where I - sqrt(-1). |
| 70 | return Double.NaN; |
| 71 | else if(z == Double.POSITIVE_INFINITY) |
| 72 | return z; |
| 73 | |
| 74 | else if(z < 0) |
| 75 | { |
| 76 | /* |
| 77 | * Using the identity |
| 78 | * |
| 79 | * __ |
| 80 | * -|| |
| 81 | * gamma(-z) = ------------------- |
| 82 | * __ |
| 83 | * z gamma(z) sin || z |
| 84 | */ |
| 85 | |
| 86 | z = -z; |
| 87 | return -PI/(z*gamma(z)*sin(PI*z)); |
| 88 | } |
| 89 | |
| 90 | /** |
| 91 | * General case for z > 0, from Numerical Recipes in C (2nd ed. Cambridge University Press, 1992). |error| is <= 2*10^-10 for all z > 0 |
| 92 | * |
| 93 | * ____ / 6 \ |
| 94 | * / __ | ===== p | |
| 95 | * \/ 2 || | \ n | z + 0.5 - (z + 5.5) |
| 96 | * gamma(z) = ------- |p + > -----| (z + 5.5) e |
| 97 | * z | 0 / z + n| |
| 98 | * | ===== | |
| 99 | * \ n = 1 / |
| 100 | * |
| 101 | * |
| 102 | * see http://www.rskey.org/gamma.htm |
| 103 | * |
| 104 | */ |
| 105 | |
| 106 | double p[] = |
| 107 | { |
| 108 | 1.000000000190015,76.18009172947146,-86.50532032941677, |
| 109 | 24.01409824083091,-1.231739572450155,1.208650973866179e-3,-5.395239384953e-6 |
| 110 | }; |
| 111 | |
| 112 | double innerLoop = 0; |
| 113 | for(int n = 1; n < p.length; n++) |
| 114 | innerLoop += p[n]/(z+n); |
| 115 | |
| 116 | double result = p[0] + innerLoop; |
| 117 | |
| 118 | result *= sqrt(2*PI)/z; |
| 119 | |
| 120 | return result*pow(z+5.5, z+0.5)*exp(-(z+5.5)); |
| 121 | } |
| 122 | |
| 123 | /** |
| 124 | * Computes the natural logarithm of {@link #gamma(double) }. |