1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 package org.apache.commons.math.special;
17
18 import java.io.Serializable;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.MathException;
22 import org.apache.commons.math.util.ContinuedFraction;
23
24 /***
25 * This is a utility class that provides computation methods related to the
26 * Gamma family of functions.
27 *
28 * @version $Revision: 1.22 $ $Date: 2004/10/08 05:53:18 $
29 */
30 public class Gamma implements Serializable {
31
32 /*** Maximum allowed numerical error. */
33 private static final double DEFAULT_EPSILON = 10e-9;
34
35 /*** Lanczos coefficients */
36 private static double[] lanczos =
37 {
38 0.99999999999999709182,
39 57.156235665862923517,
40 -59.597960355475491248,
41 14.136097974741747174,
42 -0.49191381609762019978,
43 .33994649984811888699e-4,
44 .46523628927048575665e-4,
45 -.98374475304879564677e-4,
46 .15808870322491248884e-3,
47 -.21026444172410488319e-3,
48 .21743961811521264320e-3,
49 -.16431810653676389022e-3,
50 .84418223983852743293e-4,
51 -.26190838401581408670e-4,
52 .36899182659531622704e-5,
53 };
54
55 /*** Avoid repeated computation of log of 2 PI in logGamma */
56 private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
57
58
59 /***
60 * Default constructor. Prohibit instantiation.
61 */
62 private Gamma() {
63 super();
64 }
65
66 /***
67 * Returns the natural logarithm of the gamma function Γ(x).
68 *
69 * The implementation of this method is based on:
70 * <ul>
71 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">
72 * Gamma Function</a>, equation (28).</li>
73 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
74 * Lanczos Approximation</a>, equations (1) through (5).</li>
75 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
76 * the computation of the convergent Lanczos complex Gamma approximation
77 * </a></li>
78 * </ul>
79 *
80 * @param x the value.
81 * @return log(Γ(x))
82 */
83 public static double logGamma(double x) {
84 double ret;
85
86 if (Double.isNaN(x) || (x <= 0.0)) {
87 ret = Double.NaN;
88 } else {
89 double g = 607.0 / 128.0;
90
91 double sum = 0.0;
92 for (int i = lanczos.length - 1; i > 0; --i) {
93 sum = sum + (lanczos[i] / (x + i));
94 }
95 sum = sum + lanczos[0];
96
97 double tmp = x + g + .5;
98 ret = ((x + .5) * Math.log(tmp)) - tmp +
99 HALF_LOG_2_PI + Math.log(sum / x);
100 }
101
102 return ret;
103 }
104
105 /***
106 * Returns the regularized gamma function P(a, x).
107 *
108 * @param a the a parameter.
109 * @param x the value.
110 * @return the regularized gamma function P(a, x)
111 * @throws MathException if the algorithm fails to converge.
112 */
113 public static double regularizedGammaP(double a, double x)
114 throws MathException
115 {
116 return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
117 }
118
119
120 /***
121 * Returns the regularized gamma function P(a, x).
122 *
123 * The implementation of this method is based on:
124 * <ul>
125 * <li>
126 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
127 * Regularized Gamma Function</a>, equation (1).</li>
128 * <li>
129 * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
130 * Incomplete Gamma Function</a>, equation (4).</li>
131 * <li>
132 * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
133 * Confluent Hypergeometric Function of the First Kind</a>, equation (1).
134 * </li>
135 * </ul>
136 *
137 * @param a the a parameter.
138 * @param x the value.
139 * @param epsilon When the absolute value of the nth item in the
140 * series is less than epsilon the approximation ceases
141 * to calculate further elements in the series.
142 * @param maxIterations Maximum number of "iterations" to complete.
143 * @return the regularized gamma function P(a, x)
144 * @throws MathException if the algorithm fails to converge.
145 */
146 public static double regularizedGammaP(double a,
147 double x,
148 double epsilon,
149 int maxIterations)
150 throws MathException
151 {
152 double ret;
153
154 if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
155 ret = Double.NaN;
156 } else if (x == 0.0) {
157 ret = 0.0;
158 } else if (a > 1.0 && x > a) {
159
160
161 ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
162 } else {
163
164 double n = 0.0;
165 double an = 1.0 / a;
166 double sum = an;
167 while (Math.abs(an) > epsilon && n < maxIterations) {
168
169 n = n + 1.0;
170 an = an * (x / (a + n));
171
172
173 sum = sum + an;
174 }
175 if (n >= maxIterations) {
176 throw new ConvergenceException(
177 "maximum number of iterations reached");
178 } else {
179 ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
180 }
181 }
182
183 return ret;
184 }
185
186 /***
187 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
188 *
189 * @param a the a parameter.
190 * @param x the value.
191 * @return the regularized gamma function Q(a, x)
192 * @throws MathException if the algorithm fails to converge.
193 */
194 public static double regularizedGammaQ(double a, double x)
195 throws MathException
196 {
197 return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
198 }
199
200 /***
201 * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
202 *
203 * The implementation of this method is based on:
204 * <ul>
205 * <li>
206 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
207 * Regularized Gamma Function</a>, equation (1).</li>
208 * <li>
209 * <a href=" http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
210 * Regularized incomplete gamma function: Continued fraction representations (formula 06.08.10.0003)</a></li>
211 * </ul>
212 *
213 * @param a the a parameter.
214 * @param x the value.
215 * @param epsilon When the absolute value of the nth item in the
216 * series is less than epsilon the approximation ceases
217 * to calculate further elements in the series.
218 * @param maxIterations Maximum number of "iterations" to complete.
219 * @return the regularized gamma function P(a, x)
220 * @throws MathException if the algorithm fails to converge.
221 */
222 public static double regularizedGammaQ(final double a,
223 double x,
224 double epsilon,
225 int maxIterations)
226 throws MathException
227 {
228 double ret;
229
230 if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
231 ret = Double.NaN;
232 } else if (x == 0.0) {
233 ret = 1.0;
234 } else if (x < a || a <= 1.0) {
235
236
237 ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
238 } else {
239
240 ContinuedFraction cf = new ContinuedFraction() {
241 protected double getA(int n, double x) {
242 return ((2.0 * n) + 1.0) - a + x;
243 }
244
245 protected double getB(int n, double x) {
246 return n * (a - n);
247 }
248 };
249
250 ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
251 ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * ret;
252 }
253
254 return ret;
255 }
256 }