Coverage Report - org.apache.commons.math.special.Gamma

Classes in this File	Line Coverage	Branch Coverage	Complexity
Gamma	93%	100%	2.375;2.375

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$ $Date: 2005-08-22 19:27:17 -0700 (Mon, 22 Aug 2005) $
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	32	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	32	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	0	super();
64	0	}
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	22340	if (Double.isNaN(x) || (x <= 0.0)) {
87	6	ret = Double.NaN;
88		} else {
89	22334	double g = 607.0 / 128.0;
90
91	22334	double sum = 0.0;
92	335010	for (int i = lanczos.length - 1; i > 0; --i) {
93	312676	sum = sum + (lanczos[i] / (x + i));
94		}
95	22334	sum = sum + lanczos[0];
96
97	22334	double tmp = x + g + .5;
98	22334	ret = ((x + .5) * Math.log(tmp)) - tmp +
99		HALF_LOG_2_PI + Math.log(sum / x);
100		}
101
102	22340	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	1592	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	7416	if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
155	10	ret = Double.NaN;
156	7406	} else if (x == 0.0) {
157	120	ret = 0.0;
158	7286	} else if (a >= 1.0 && x > a) {
159		// use regularizedGammaQ because it should converge faster in this
160		// case.
161	604	ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
162		} else {
163		// calculate series
164	6682	double n = 0.0; // current element index
165	6682	double an = 1.0 / a; // n-th element in the series
166	6682	double sum = an; // partial sum
167	5618004	while (Math.abs(an) > epsilon && n < maxIterations) {
168		// compute next element in the series
169	5611322	n = n + 1.0;
170	5611322	an = an * (x / (a + n));
171
172		// update partial sum
173	5611322	sum = sum + an;
174		}
175	6682	if (n >= maxIterations) {
176	0	throw new ConvergenceException(
177		"maximum number of iterations reached");
178		} else {
179	6682	ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
180		}
181		}
182
183	7416	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	14	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	6964	if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
231	10	ret = Double.NaN;
232	6954	} else if (x == 0.0) {
233	2	ret = 1.0;
234	6952	} else if (x < a || a < 1.0) {
235		// use regularizedGammaP because it should converge faster in this
236		// case.
237	5354	ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
238		} else {
239		// create continued fraction
240	1598	ContinuedFraction cf = new ContinuedFraction() {
241		protected double getA(int n, double x) {
242	39084	return ((2.0 * n) + 1.0) - a + x;
243		}
244
245	1598	protected double getB(int n, double x) {
246	37486	return n * (a - n);
247		}
248		};
249
250	1598	ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
251	1598	ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * ret;
252		}
253
254	6964	return ret;
255		}
256		}