001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.numbers.gamma; 018 019import java.text.MessageFormat; 020 021import org.apache.commons.numbers.fraction.ContinuedFraction; 022 023/** 024 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 025 * Regularized Gamma functions</a>. 026 * 027 * Class is immutable. 028 */ 029public final class RegularizedGamma { 030 /** Maximum allowed numerical error. */ 031 private static final double DEFAULT_EPSILON = 1e-15; 032 033 /** Private constructor. */ 034 private RegularizedGamma() { 035 // intentionally empty. 036 } 037 038 /** 039 * \( P(a, x) \) <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 040 * regularized Gamma function</a>. 041 * 042 * Class is immutable. 043 */ 044 public static final class P { 045 /** Prevent instantiation. */ 046 private P() {} 047 048 /** 049 * Computes the regularized gamma function \( P(a, x) \). 050 * 051 * @param a Argument. 052 * @param x Argument. 053 * @return \( P(a, x) \). 054 * @throws ArithmeticException if the continued fraction fails to converge. 055 */ 056 public static double value(double a, 057 double x) { 058 return value(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); 059 } 060 061 /** 062 * Computes the regularized gamma function \( P(a, x) \). 063 * 064 * The implementation of this method is based on: 065 * <ul> 066 * <li> 067 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 068 * Regularized Gamma Function</a>, equation (1) 069 * </li> 070 * <li> 071 * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html"> 072 * Incomplete Gamma Function</a>, equation (4). 073 * </li> 074 * <li> 075 * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html"> 076 * Confluent Hypergeometric Function of the First Kind</a>, equation (1). 077 * </li> 078 * </ul> 079 * 080 * @param a Argument. 081 * @param x Argument. 082 * @param epsilon Tolerance in continued fraction evaluation. 083 * @param maxIterations Maximum number of iterations in continued fraction evaluation. 084 * @return \( P(a, x) \). 085 * @throws ArithmeticException if the continued fraction fails to converge. 086 */ 087 public static double value(double a, 088 double x, 089 double epsilon, 090 int maxIterations) { 091 if (Double.isNaN(a) || 092 Double.isNaN(x) || 093 a <= 0 || 094 x < 0) { 095 return Double.NaN; 096 } else if (x == 0) { 097 return 0; 098 } else if (x >= a + 1) { 099 // Q should converge faster in this case. 100 return 1 - RegularizedGamma.Q.value(a, x, epsilon, maxIterations); 101 } else { 102 // Series. 103 double n = 0; // current element index 104 double an = 1 / a; // n-th element in the series 105 double sum = an; // partial sum 106 while (Math.abs(an / sum) > epsilon && 107 n < maxIterations && 108 sum < Double.POSITIVE_INFINITY) { 109 // compute next element in the series 110 n += 1; 111 an *= x / (a + n); 112 113 // update partial sum 114 sum += an; 115 } 116 if (n >= maxIterations) { 117 throw new ArithmeticException( 118 MessageFormat.format("Failed to converge within {0} iterations", maxIterations)); 119 } else if (Double.isInfinite(sum)) { 120 return 1; 121 } else { 122 // Ensure result is in the range [0, 1] 123 final double result = Math.exp(-x + (a * Math.log(x)) - LogGamma.value(a)) * sum; 124 return result > 1.0 ? 1.0 : result; 125 } 126 } 127 } 128 } 129 130 /** 131 * Creates the \( Q(a, x) \equiv 1 - P(a, x) \) <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 132 * regularized Gamma function</a>. 133 * 134 * Class is immutable. 135 */ 136 public static final class Q { 137 /** Prevent instantiation. */ 138 private Q() {} 139 140 /** 141 * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \). 142 * 143 * @param a Argument. 144 * @param x Argument. 145 * @return \( Q(a, x) \). 146 * @throws ArithmeticException if the continued fraction fails to converge. 147 */ 148 public static double value(double a, 149 double x) { 150 return value(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); 151 } 152 153 /** 154 * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \). 155 * 156 * The implementation of this method is based on: 157 * <ul> 158 * <li> 159 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html"> 160 * Regularized Gamma Function</a>, equation (1). 161 * </li> 162 * <li> 163 * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/"> 164 * Regularized incomplete gamma function: Continued fraction representations 165 * (formula 06.08.10.0003)</a> 166 * </li> 167 * </ul> 168 * 169 * @param a Argument. 170 * @param x Argument. 171 * @param epsilon Tolerance in continued fraction evaluation. 172 * @param maxIterations Maximum number of iterations in continued fraction evaluation. 173 * @throws ArithmeticException if the continued fraction fails to converge. 174 * @return \( Q(a, x) \). 175 */ 176 public static double value(final double a, 177 double x, 178 double epsilon, 179 int maxIterations) { 180 if (Double.isNaN(a) || 181 Double.isNaN(x) || 182 a <= 0 || 183 x < 0) { 184 return Double.NaN; 185 } else if (x == 0) { 186 return 1; 187 } else if (x < a + 1) { 188 // P should converge faster in this case. 189 return 1 - RegularizedGamma.P.value(a, x, epsilon, maxIterations); 190 } else { 191 final ContinuedFraction cf = new ContinuedFraction() { 192 /** {@inheritDoc} */ 193 @Override 194 protected double getA(int n, double x) { 195 return n * (a - n); 196 } 197 198 /** {@inheritDoc} */ 199 @Override 200 protected double getB(int n, double x) { 201 return ((2 * n) + 1) - a + x; 202 } 203 }; 204 205 return Math.exp(-x + (a * Math.log(x)) - LogGamma.value(a)) / 206 cf.evaluate(x, epsilon, maxIterations); 207 } 208 } 209 } 210}