View Javadoc

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;
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(&#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 }