KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > org > apache > commons > math > special > Gamma


1 /*
2  * Copyright 2003-2004 The Apache Software Foundation.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 package org.apache.commons.math.special;
17
18 import java.io.Serializable JavaDoc;
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$ $Date: 2005-02-26 05:11:52 -0800 (Sat, 26 Feb 2005) $
29  */
30 public class Gamma implements Serializable JavaDoc {
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(&#915;(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             // use regularizedGammaQ because it should converge faster in this
160 // case.
161 ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
162         } else {
163             // calculate series
164 double n = 0.0; // current element index
165 double an = 1.0 / a; // n-th element in the series
166 double sum = an; // partial sum
167 while (Math.abs(an) > epsilon && n < maxIterations) {
168                 // compute next element in the series
169 n = n + 1.0;
170                 an = an * (x / (a + n));
171
172                 // update partial sum
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             // use regularizedGammaP because it should converge faster in this
236 // case.
237 ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
238         } else {
239             // create continued fraction
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 }
257
Popular Tags