Java tutorial
/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.special; import org.apache.commons.math.MathException; import org.apache.commons.math.MaxIterationsExceededException; import org.apache.commons.math.util.ContinuedFraction; import org.apache.commons.math.util.FastMath; /** * This is a utility class that provides computation methods related to the * Gamma family of functions. * * @version $Revision: 1042510 $ $Date: 2010-12-06 02:54:18 +0100 (lun. 06 dc. 2010) $ */ public class Gamma { /** * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a> * @since 2.0 */ public static final double GAMMA = 0.577215664901532860606512090082; /** Maximum allowed numerical error. */ private static final double DEFAULT_EPSILON = 10e-15; /** Lanczos coefficients */ private static final double[] LANCZOS = { 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248, 14.136097974741747174, -0.49191381609762019978, .33994649984811888699e-4, .46523628927048575665e-4, -.98374475304879564677e-4, .15808870322491248884e-3, -.21026444172410488319e-3, .21743961811521264320e-3, -.16431810653676389022e-3, .84418223983852743293e-4, -.26190838401581408670e-4, .36899182659531622704e-5, }; /** Avoid repeated computation of log of 2 PI in logGamma */ private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI); // limits for switching algorithm in digamma /** C limit. */ private static final double C_LIMIT = 49; /** S limit. */ private static final double S_LIMIT = 1e-5; /** * Default constructor. Prohibit instantiation. */ private Gamma() { super(); } /** * Returns the natural logarithm of the gamma function Γ(x). * * The implementation of this method is based on: * <ul> * <li><a href="http://mathworld.wolfram.com/GammaFunction.html"> * Gamma Function</a>, equation (28).</li> * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html"> * Lanczos Approximation</a>, equations (1) through (5).</li> * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on * the computation of the convergent Lanczos complex Gamma approximation * </a></li> * </ul> * * @param x the value. * @return log(Γ(x)) */ public static double logGamma(double x) { double ret; if (Double.isNaN(x) || (x <= 0.0)) { ret = Double.NaN; } else { double g = 607.0 / 128.0; double sum = 0.0; for (int i = LANCZOS.length - 1; i > 0; --i) { sum = sum + (LANCZOS[i] / (x + i)); } sum = sum + LANCZOS[0]; double tmp = x + g + .5; ret = ((x + .5) * FastMath.log(tmp)) - tmp + HALF_LOG_2_PI + FastMath.log(sum / x); } return ret; } /** * Returns the regularized gamma function P(a, x). * * @param a the a parameter. * @param x the value. * @return the regularized gamma function P(a, x) * @throws MathException if the algorithm fails to converge. */ public static double regularizedGammaP(double a, double x) throws MathException { return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); } /** * Returns the regularized gamma function P(a, x). * * The implementation of this method is based on: * <ul> * <li> * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> * Regularized Gamma Function</a>, equation (1).</li> * <li> * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html"> * Incomplete Gamma Function</a>, equation (4).</li> * <li> * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html"> * Confluent Hypergeometric Function of the First Kind</a>, equation (1). * </li> * </ul> * * @param a the a parameter. * @param x the value. * @param epsilon When the absolute value of the nth item in the * series is less than epsilon the approximation ceases * to calculate further elements in the series. * @param maxIterations Maximum number of "iterations" to complete. * @return the regularized gamma function P(a, x) * @throws MathException if the algorithm fails to converge. */ public static double regularizedGammaP(double a, double x, double epsilon, int maxIterations) throws MathException { double ret; if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { ret = Double.NaN; } else if (x == 0.0) { ret = 0.0; } else if (x >= a + 1) { // use regularizedGammaQ because it should converge faster in this // case. ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations); } else { // calculate series double n = 0.0; // current element index double an = 1.0 / a; // n-th element in the series double sum = an; // partial sum while (FastMath.abs(an / sum) > epsilon && n < maxIterations && sum < Double.POSITIVE_INFINITY) { // compute next element in the series n = n + 1.0; an = an * (x / (a + n)); // update partial sum sum = sum + an; } if (n >= maxIterations) { throw new MaxIterationsExceededException(maxIterations); } else if (Double.isInfinite(sum)) { ret = 1.0; } else { ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum; } } return ret; } /** * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). * * @param a the a parameter. * @param x the value. * @return the regularized gamma function Q(a, x) * @throws MathException if the algorithm fails to converge. */ public static double regularizedGammaQ(double a, double x) throws MathException { return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); } /** * Returns the regularized gamma function Q(a, x) = 1 - P(a, x). * * The implementation of this method is based on: * <ul> * <li> * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> * Regularized Gamma Function</a>, equation (1).</li> * <li> * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/"> * Regularized incomplete gamma function: Continued fraction representations (formula 06.08.10.0003)</a></li> * </ul> * * @param a the a parameter. * @param x the value. * @param epsilon When the absolute value of the nth item in the * series is less than epsilon the approximation ceases * to calculate further elements in the series. * @param maxIterations Maximum number of "iterations" to complete. * @return the regularized gamma function P(a, x) * @throws MathException if the algorithm fails to converge. */ public static double regularizedGammaQ(final double a, double x, double epsilon, int maxIterations) throws MathException { double ret; if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { ret = Double.NaN; } else if (x == 0.0) { ret = 1.0; } else if (x < a + 1.0) { // use regularizedGammaP because it should converge faster in this // case. ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations); } else { // create continued fraction ContinuedFraction cf = new ContinuedFraction() { @Override protected double getA(int n, double x) { return ((2.0 * n) + 1.0) - a + x; } @Override protected double getB(int n, double x) { return n * (a - n); } }; ret = 1.0 / cf.evaluate(x, epsilon, maxIterations); ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret; } return ret; } /** * <p>Computes the digamma function of x.</p> * * <p>This is an independently written implementation of the algorithm described in * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p> * * <p>Some of the constants have been changed to increase accuracy at the moderate expense * of run-time. The result should be accurate to within 10^-8 absolute tolerance for * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p> * * <p>Performance for large negative values of x will be quite expensive (proportional to * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results * less than 10^5 and 10^-8 relative for results larger than that.</p> * * @param x the argument * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller * @see <a href="http://en.wikipedia.org/wiki/Digamma_function"> Digamma at wikipedia </a> * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf"> Bernardo's original article </a> * @since 2.0 */ public static double digamma(double x) { if (x > 0 && x <= S_LIMIT) { // use method 5 from Bernardo AS103 // accurate to O(x) return -GAMMA - 1 / x; } if (x >= C_LIMIT) { // use method 4 (accurate to O(1/x^8) double inv = 1 / (x * x); // 1 1 1 1 // log(x) - --- - ------ + ------- - ------- // 2 x 12 x^2 120 x^4 252 x^6 return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); } return digamma(x + 1) - 1 / x; } /** * <p>Computes the trigamma function of x. This function is derived by taking the derivative of * the implementation of digamma.</p> * * @param x the argument * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function"> Trigamma at wikipedia </a> * @see Gamma#digamma(double) * @since 2.0 */ public static double trigamma(double x) { if (x > 0 && x <= S_LIMIT) { return 1 / (x * x); } if (x >= C_LIMIT) { double inv = 1 / (x * x); // 1 1 1 1 1 // - + ---- + ---- - ----- + ----- // x 2 3 5 7 // 2 x 6 x 30 x 42 x return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); } return trigamma(x + 1) + 1 / (x * x); } }