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