Class BoostGamma
This code has been adapted from the Boost
c++
implementation <boost/math/special_functions/gamma.hpp>
.
All work is copyright to the original authors and subject to the Boost Software License.
- See Also:
-
Nested Class Summary
Nested ClassesModifier and TypeClassDescription(package private) static final class
53-bit precision implementation of the Lanczos approximation. -
Field Summary
FieldsModifier and TypeFieldDescriptionprivate static final double
Default epsilon value for relative error.private static final double
Euler's constant.private static final double[]
All factorials that can be represented as a double.private static final int
The threshold value for choosing the Lanczos approximation.private static final int
Approximate value for ln(Double.MAX_VALUE).private static final int
Approximate value for ln(Double.MIN_VALUE).private static final double
ln(pi).private static final double
ln(sqrt(2 pi)).private static final int
The largest factorial that can be represented as a double.private static final int
The largest integer value for gamma(z) that can be represented as a double.private static final double
Value for the sqrt of the epsilon for relative error.private static final double
2^53. -
Constructor Summary
Constructors -
Method Summary
Modifier and TypeMethodDescriptionprivate static double
finiteGammaQ
(double a, double x) Upper gamma fraction for integer a.private static double
finiteHalfGammaQ
(double a, double x) Upper gamma fraction for half integer a.(package private) static double
fullIgammaPrefix
(double a, double z) Calculate power term prefix (z^a)(e^-z) used in the non-normalised incomplete gammas.private static double
gammaIncompleteImp
(double a, double x, boolean normalised, boolean invert, Policy pol) Main incomplete gamma entry point, handles all four incomplete gammas.(package private) static double
gammaP
(double a, double x) Regularised lower incomplete gamma.(package private) static double
Regularised lower incomplete gamma.(package private) static double
gammaPDerivative
(double a, double x) Derivative of the regularised lower incomplete gamma.(package private) static double
gammaQ
(double a, double x) Regularised upper incomplete gamma.(package private) static double
Regularised upper incomplete gamma.(package private) static double[]
All factorials that are representable as a double.(package private) static double
igammaTemmeLarge
(double a, double x) Implements the asymptotic expansions of the incomplete gamma functions P(a, x) and Q(a, x), used when a is large and x ~ a.(package private) static double
incompleteTgammaLargeX
(double a, double x, Policy pol) Incomplete tgamma for large X.private static boolean
isOdd
(double v) Checks if the value is odd.(package private) static double
lgamma
(double z) Log Gamma function.(package private) static double
lgamma
(double z, int[] sign) Log Gamma function.private static double
lgammaSmall
(double z, double zm1, double zm2) Log Gamma function for small z.(package private) static double
lowerGammaSeries
(double a, double z, double initValue, Policy pol) Lower gamma series.private static boolean
nonZeroLength
(int[] array) Return true if the array is not null and has non-zero length.(package private) static double
regularisedGammaPrefix
(double a, double z) Compute (z^a)(e^-z)/tgamma(a).(package private) static double
sinpx
(double x) Ad hoc function calculates x * sin(pi * x), taking extra care near when x is near a whole number.(package private) static double
tgamma
(double z) Gamma function.(package private) static double
tgamma
(double a, double x) Full upper incomplete gamma.(package private) static double
Full upper incomplete gamma.(package private) static double
tgamma1pm1
(double dz) Calculates tgamma(1+dz)-1.(package private) static double
tgammaDeltaRatio
(double z, double delta) Ratio of gamma functions.private static double
tgammaDeltaRatioImpLanczos
(double z, double delta) Ratio of gamma functions using Lanczos support.(package private) static double
tgammaLower
(double a, double x) Full lower incomplete gamma.(package private) static double
tgammaLower
(double a, double x, Policy policy) Full lower incomplete gamma.(package private) static double
tgammaRatio
(double x, double y) Ratio of gamma functions.private static double
tgammaSmallUpperPart
(double a, double x, Policy pol, double[] pgam, boolean invert) Upper gamma fraction for very small a.(package private) static double
uncheckedFactorial
(int n) Returns the factorial of n.(package private) static double
upperGammaFraction
(double a, double z, Policy pol) Upper gamma fraction.
-
Field Details
-
EPSILON
private static final double EPSILONDefault epsilon value for relative error. This is equal to the Boost constantboost::math::tools::epsilon<double>()
.- See Also:
-
ROOT_EPSILON
private static final double ROOT_EPSILONValue for the sqrt of the epsilon for relative error. This is equal to the Boost constantboost::math::tools::root_epsilon<double>()
.- See Also:
-
LOG_MAX_VALUE
private static final int LOG_MAX_VALUEApproximate value for ln(Double.MAX_VALUE). This is equal to the Boost constantboost::math::tools::log_max_value<double>()
. No termx
should be used inexp(x)
ifx > LOG_MAX_VALUE
to avoid overflow.- See Also:
-
LOG_MIN_VALUE
private static final int LOG_MIN_VALUEApproximate value for ln(Double.MIN_VALUE). This is equal to the Boost constantboost::math::tools::log_min_value<double>()
. No termx
should be used inexp(x)
ifx < LOG_MIN_VALUE
to avoid underflow to sub-normal or zero.- See Also:
-
MAX_FACTORIAL
private static final int MAX_FACTORIALThe largest factorial that can be represented as a double. This is equal to the Boost constantboost::math::max_factorial<double>::value
.- See Also:
-
MAX_GAMMA_Z
private static final int MAX_GAMMA_ZThe largest integer value for gamma(z) that can be represented as a double.- See Also:
-
LOG_ROOT_TWO_PI
private static final double LOG_ROOT_TWO_PIln(sqrt(2 pi)). Computed to 25-digits precision.- See Also:
-
LOG_PI
private static final double LOG_PIln(pi). Computed to 25-digits precision.- See Also:
-
EULER
private static final double EULEREuler's constant.- See Also:
-
LANCZOS_THRESHOLD
private static final int LANCZOS_THRESHOLDThe threshold value for choosing the Lanczos approximation.- See Also:
-
TWO_POW_53
private static final double TWO_POW_532^53.- See Also:
-
FACTORIAL
private static final double[] FACTORIALAll factorials that can be represented as a double. Size = 171.
-
-
Constructor Details
-
BoostGamma
private BoostGamma()Private constructor.
-
-
Method Details
-
getFactorials
static double[] getFactorials()All factorials that are representable as a double. This data is exposed for testing.- Returns:
- factorials
-
uncheckedFactorial
static double uncheckedFactorial(int n) Returns the factorial of n. This is unchecked as an index out of bound exception will occur if the value n is not within the range [0, 170]. This function is exposed for use inBoostBeta
.- Parameters:
n
- Argument n (must be in [0, 170])- Returns:
- n!
-
tgamma
static double tgamma(double z) Gamma function.For small
z
this is based on the NSWC Library of Mathematics Subroutines double precision implementation,DGAMMA
.For large
z
this is an implementation of the Boost C++ tgamma function with Lanczos support.Integers are handled using a look-up table of factorials.
Note: The Boost C++ implementation uses the Lanczos sum for all
z
. When promotion of double to long double is not available this has larger errors than the double precision specific NSWC implementation. For largerz
the Boost C++ Lanczos implementation incorporates the sqrt(2 pi) factor and has lower error than the implementation using theLanczosApproximation
class.- Parameters:
z
- Argument z- Returns:
- gamma value
-
sinpx
static double sinpx(double x) Ad hoc function calculates x * sin(pi * x), taking extra care near when x is near a whole number.- Parameters:
x
- Value (assumed to be negative)- Returns:
- x * sin(pi * x)
-
isOdd
private static boolean isOdd(double v) Checks if the value is odd.- Parameters:
v
- Value (assumed to be positive and an integer)- Returns:
- true if odd
-
lgamma
static double lgamma(double z) Log Gamma function. Defined as the natural logarithm of the absolute value of tgamma(z).- Parameters:
z
- Argument z- Returns:
- log gamma value
-
lgamma
static double lgamma(double z, int[] sign) Log Gamma function. Defined as the natural logarithm of the absolute value of tgamma(z).- Parameters:
z
- Argument zsign
- If a non-zero length array the first index is set on output to the sign of tgamma(z)- Returns:
- log gamma value
-
lgammaSmall
private static double lgammaSmall(double z, double zm1, double zm2) Log Gamma function for small z.- Parameters:
z
- Argument zzm1
-z - 1
zm2
-z - 2
- Returns:
- log gamma value
-
tgamma1pm1
static double tgamma1pm1(double dz) Calculates tgamma(1+dz)-1.- Parameters:
dz
- Argument- Returns:
- tgamma(1+dz)-1
-
tgamma
static double tgamma(double a, double x) Full upper incomplete gamma.- Parameters:
a
- Argument ax
- Argument x- Returns:
- upper gamma value
-
tgamma
Full upper incomplete gamma.- Parameters:
a
- Argument ax
- Argument xpolicy
- Function evaluation policy- Returns:
- upper gamma value
-
tgammaLower
static double tgammaLower(double a, double x) Full lower incomplete gamma.- Parameters:
a
- Argument ax
- Argument x- Returns:
- lower gamma value
-
tgammaLower
Full lower incomplete gamma.- Parameters:
a
- Argument ax
- Argument xpolicy
- Function evaluation policy- Returns:
- lower gamma value
-
gammaQ
static double gammaQ(double a, double x) Regularised upper incomplete gamma.- Parameters:
a
- Argument ax
- Argument x- Returns:
- q
-
gammaQ
Regularised upper incomplete gamma.- Parameters:
a
- Argument ax
- Argument xpolicy
- Function evaluation policy- Returns:
- q
-
gammaP
static double gammaP(double a, double x) Regularised lower incomplete gamma.- Parameters:
a
- Argument ax
- Argument x- Returns:
- p
-
gammaP
Regularised lower incomplete gamma.- Parameters:
a
- Argument ax
- Argument xpolicy
- Function evaluation policy- Returns:
- p
-
gammaPDerivative
static double gammaPDerivative(double a, double x) Derivative of the regularised lower incomplete gamma.\( \frac{e^{-x} x^{a-1}}{\Gamma(a)} \)
Adapted from
boost::math::detail::gamma_p_derivative_imp
- Parameters:
a
- Argument ax
- Argument x- Returns:
- p derivative
-
gammaIncompleteImp
private static double gammaIncompleteImp(double a, double x, boolean normalised, boolean invert, Policy pol) Main incomplete gamma entry point, handles all four incomplete gammas. Adapted fromboost::math::detail::gamma_incomplete_imp
.The Boost code has a pointer
p_derivative
that can be set to the value of the derivative. This is used for the inverse incomplete gamma functionsgamma_(p|q)_inv_imp
. It is not required for the forward evaluation functions.- Parameters:
a
- Argument ax
- Argument xnormalised
- true to compute the regularised valueinvert
- true to compute the upper value Q (default is lower value P)pol
- Function evaluation policy- Returns:
- gamma value
-
upperGammaFraction
Upper gamma fraction. Multiply result by z^a * e^-z to get the full upper incomplete integral. Divide by tgamma(z) to normalise.- Parameters:
a
- Argument az
- Argument zpol
- Function evaluation policy- Returns:
- upper gamma fraction
-
finiteGammaQ
private static double finiteGammaQ(double a, double x) Upper gamma fraction for integer a. Called whena < 30
and-x > LOG_MIN_VALUE
.- Parameters:
a
- Argument a (assumed to be small)x
- Argument x- Returns:
- upper gamma fraction
-
finiteHalfGammaQ
private static double finiteHalfGammaQ(double a, double x) Upper gamma fraction for half integer a. Called whena < 30
and-x > LOG_MIN_VALUE
.- Parameters:
a
- Argument a (assumed to be small)x
- Argument x- Returns:
- upper gamma fraction
-
lowerGammaSeries
Lower gamma series. Multiply result by ((z^a) * (e^-z) / a) to get the full lower incomplete integral. Then divide by tgamma(a) to get the normalised value.- Parameters:
a
- Argument az
- Argument zinitValue
- Initial valuepol
- Function evaluation policy- Returns:
- lower gamma series
-
tgammaSmallUpperPart
private static double tgammaSmallUpperPart(double a, double x, Policy pol, double[] pgam, boolean invert) Upper gamma fraction for very small a.- Parameters:
a
- Argument a (assumed to be small)x
- Argument xpol
- Function evaluation policypgam
- set to value of gamma(a) on outputinvert
- true to invert the result- Returns:
- upper gamma fraction
-
fullIgammaPrefix
static double fullIgammaPrefix(double a, double z) Calculate power term prefix (z^a)(e^-z) used in the non-normalised incomplete gammas.- Parameters:
a
- Argument az
- Argument z- Returns:
- incomplete gamma prefix
-
regularisedGammaPrefix
static double regularisedGammaPrefix(double a, double z) Compute (z^a)(e^-z)/tgamma(a).Most of the error occurs in this function.
- Parameters:
a
- Argument az
- Argument z- Returns:
- regularized gamma prefix
-
igammaTemmeLarge
static double igammaTemmeLarge(double a, double x) Implements the asymptotic expansions of the incomplete gamma functions P(a, x) and Q(a, x), used when a is large and x ~ a.The primary reference is:
"The Asymptotic Expansion of the Incomplete Gamma Functions" N. M. Temme. Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
A different way of evaluating these expansions, plus a lot of very useful background information is in:
"A Set of Algorithms For the Incomplete Gamma Functions." N. M. Temme. Probability in the Engineering and Informational Sciences, 8, 1994, 291.
An alternative implementation is in:
"Computation of the Incomplete Gamma Function Ratios and their Inverse." A. R. Didonato and A. H. Morris. ACM TOMS, Vol 12, No 4, Dec 1986, p377.
This is a port of the function accurate for 53-bit mantissas (IEEE double precision or 10^-17). To understand the code, refer to Didonato and Morris, from Eq 17 and 18 onwards.
The coefficients used here are not taken from Didonato and Morris: the domain over which these expansions are used is slightly different to theirs, and their constants are not quite accurate enough for 128-bit long doubles. Instead the coefficients were calculated using the methods described by Temme p762 from Eq 3.8 onwards. The values obtained agree with those obtained by Didonato and Morris (at least to the first 30 digits that they provide). At double precision the degrees of polynomial required for full machine precision are close to those recommended to Didonato and Morris, but of course many more terms are needed for larger types.
Adapted from
boost/math/special_functions/detail/igamma_large.hpp
.- Parameters:
a
- the ax
- the x- Returns:
- the double
-
incompleteTgammaLargeX
Incomplete tgamma for large X.This summation is a source of error as the series starts at 1 and descends to zero. It can have thousands of iterations when a and z are large and close in value.
- Parameters:
a
- Argument ax
- Argument xpol
- Function evaluation policy- Returns:
- incomplete tgamma
-
nonZeroLength
private static boolean nonZeroLength(int[] array) Return true if the array is not null and has non-zero length.- Parameters:
array
- Array- Returns:
- true if a non-zero length array
-
tgammaDeltaRatio
static double tgammaDeltaRatio(double z, double delta) Ratio of gamma functions.\[ tgamma_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)} \]
Adapted from
tgamma_delta_ratio_imp
. The use ofmax_factorial<double>::value == 170
has been replaced withMAX_GAMMA_Z == 171
. This threshold is used when it is possible to call the gamma function without overflow.- Parameters:
z
- Argument zdelta
- The difference- Returns:
- gamma ratio
-
tgammaDeltaRatioImpLanczos
private static double tgammaDeltaRatioImpLanczos(double z, double delta) Ratio of gamma functions using Lanczos support.\[ tgamma_delta_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)}
Adapted from
tgamma_delta_ratio_imp_lanczos
. The use ofmax_factorial<double>::value == 170
has been replaced withMAX_GAMMA_Z == 171
. This threshold is used when it is possible to use the precomputed factorial table.- Parameters:
z
- Argument zdelta
- The difference- Returns:
- gamma ratio
-
tgammaRatio
static double tgammaRatio(double x, double y) Ratio of gamma functions.\[ tgamma_ratio(x, y) = \frac{\Gamma(x)}{\Gamma(y)}
Adapted from
tgamma_ratio_imp
. The use ofmax_factorial<double>::value == 170
has been replaced withMAX_GAMMA_Z == 171
. This threshold is used when it is possible to call the gamma function without overflow.- Parameters:
x
- Argument x (must be positive finite)y
- Argument y (must be positive finite)- Returns:
- gamma ratio (or nan)
-