in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java [1185:1466]
private static double gammaIncompleteImp(double a, double x,
boolean normalised, boolean invert, Policy pol) {
if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
return Double.NaN;
}
double result = 0;
if (a >= MAX_FACTORIAL && !normalised) {
//
// When we're computing the non-normalized incomplete gamma
// and a is large the result is rather hard to compute unless
// we use logs. There are really two options - if x is a long
// way from a in value then we can reliably use methods 2 and 4
// below in logarithmic form and go straight to the result.
// Otherwise we let the regularized gamma take the strain
// (the result is unlikely to underflow in the central region anyway)
// and combine with lgamma in the hopes that we get a finite result.
//
if (invert && (a * 4 < x)) {
// This is method 4 below, done in logs:
result = a * Math.log(x) - x;
result += Math.log(upperGammaFraction(a, x, pol));
} else if (!invert && (a > 4 * x)) {
// This is method 2 below, done in logs:
result = a * Math.log(x) - x;
result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
} else {
result = gammaIncompleteImp(a, x, true, invert, pol);
if (result == 0) {
if (invert) {
// Try http://functions.wolfram.com/06.06.06.0039.01
result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
result = Math.log(result) - a + (a - 0.5f) * Math.log(a) + LOG_ROOT_TWO_PI;
} else {
// This is method 2 below, done in logs, we're really outside the
// range of this method, but since the result is almost certainly
// infinite, we should probably be OK:
result = a * Math.log(x) - x;
result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
}
} else {
result = Math.log(result) + lgamma(a);
}
}
// If result is > log(MAX_VALUE) the result will overflow.
// There is no exception for this case so just return the result.
return Math.exp(result);
}
final boolean isInt;
final boolean isHalfInt;
// Update. x must be safe for exp(-x). Change to -x > LOG_MIN_VALUE.
final boolean isSmallA = (a < 30) && (a <= x + 1) && (-x > LOG_MIN_VALUE);
if (isSmallA) {
final double fa = Math.floor(a);
isInt = fa == a;
isHalfInt = !isInt && (Math.abs(fa - a) == 0.5f);
} else {
isInt = isHalfInt = false;
}
final int evalMethod;
if (isInt && (x > 0.6)) {
// calculate Q via finite sum:
invert = !invert;
evalMethod = 0;
} else if (isHalfInt && (x > 0.2)) {
// calculate Q via finite sum for half integer a:
invert = !invert;
evalMethod = 1;
} else if ((x < ROOT_EPSILON) && (a > 1)) {
evalMethod = 6;
} else if ((x > 1000) && (a < x * 0.75f)) {
// Note:
// The branch is used in Boost when:
// ((x > 1000) && ((a < x) || (Math.abs(a - 50) / x < 1)))
//
// This case was added after Boost 1_68_0.
// See: https://github.com/boostorg/math/issues/168
//
// When using only double precision for the evaluation
// it is a source of error when a ~ z as the asymptotic approximation
// sums terms t_n+1 = t_n * (a - n - 1) / z starting from t_0 = 1.
// These terms are close to 1 when a ~ z and the sum has many terms
// with reduced precision.
// This has been updated to allow only cases with fast convergence.
// It will be used when x -> infinity and a << x.
// calculate Q via asymptotic approximation:
invert = !invert;
evalMethod = 7;
} else if (x < 0.5) {
//
// Changeover criterion chosen to give a changeover at Q ~ 0.33
//
if (-0.4 / Math.log(x) < a) {
// Compute P
evalMethod = 2;
} else {
evalMethod = 3;
}
} else if (x < 1.1) {
//
// Changover here occurs when P ~ 0.75 or Q ~ 0.25:
//
if (x * 0.75f < a) {
// Compute P
evalMethod = 2;
} else {
evalMethod = 3;
}
} else {
//
// Begin by testing whether we're in the "bad" zone
// where the result will be near 0.5 and the usual
// series and continued fractions are slow to converge:
//
boolean useTemme = false;
if (normalised && (a > 20)) {
final double sigma = Math.abs((x - a) / a);
if (a > 200) {
//
// This limit is chosen so that we use Temme's expansion
// only if the result would be larger than about 10^-6.
// Below that the regular series and continued fractions
// converge OK, and if we use Temme's method we get increasing
// errors from the dominant erfc term as its (inexact) argument
// increases in magnitude.
//
if (20 / a > sigma * sigma) {
useTemme = true;
}
} else {
// Note in this zone we can't use Temme's expansion for
// types longer than an 80-bit real:
// it would require too many terms in the polynomials.
if (sigma < 0.4) {
useTemme = true;
}
}
}
if (useTemme) {
evalMethod = 5;
} else {
//
// Regular case where the result will not be too close to 0.5.
//
// Changeover here occurs at P ~ Q ~ 0.5
// Note that series computation of P is about x2 faster than continued fraction
// calculation of Q, so try and use the CF only when really necessary,
// especially for small x.
//
if (x - (1 / (3 * x)) < a) {
evalMethod = 2;
} else {
evalMethod = 4;
invert = !invert;
}
}
}
switch (evalMethod) {
case 0:
result = finiteGammaQ(a, x);
if (!normalised) {
result *= tgamma(a);
}
break;
case 1:
result = finiteHalfGammaQ(a, x);
if (!normalised) {
result *= tgamma(a);
}
break;
case 2:
// Compute P:
result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
if (result != 0) {
//
// If we're going to be inverting the result then we can
// reduce the number of series evaluations by quite
// a few iterations if we set an initial value for the
// series sum based on what we'll end up subtracting it from
// at the end.
// Have to be careful though that this optimization doesn't
// lead to spurious numeric overflow. Note that the
// scary/expensive overflow checks below are more often
// than not bypassed in practice for "sensible" input
// values:
//
double initValue = 0;
boolean optimisedInvert = false;
if (invert) {
initValue = normalised ? 1 : tgamma(a);
if (normalised || (result >= 1) || (Double.MAX_VALUE * result > initValue)) {
initValue /= result;
if (normalised || (a < 1) || (Double.MAX_VALUE / a > initValue)) {
initValue *= -a;
optimisedInvert = true;
} else {
initValue = 0;
}
} else {
initValue = 0;
}
}
result *= lowerGammaSeries(a, x, initValue, pol) / a;
if (optimisedInvert) {
invert = false;
result = -result;
}
}
break;
case 3:
// Compute Q:
invert = !invert;
final double[] g = {0};
result = tgammaSmallUpperPart(a, x, pol, g, invert);
invert = false;
if (normalised) {
// Addition to the Boost code:
if (g[0] == Double.POSITIVE_INFINITY) {
// Very small a will overflow gamma(a). Resort to logs.
// This method requires improvement as the error is very large.
// It is better than returning zero for a non-zero result.
result = Math.exp(Math.log(result) - lgamma(a));
} else {
result /= g[0];
}
}
break;
case 4:
// Compute Q:
result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
if (result != 0) {
result *= upperGammaFraction(a, x, pol);
}
break;
case 5:
// Call 53-bit version
result = igammaTemmeLarge(a, x);
if (x >= a) {
invert = !invert;
}
break;
case 6:
// x is so small that P is necessarily very small too, use
// http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
if (normalised) {
// If tgamma overflows then result = 0
result = Math.pow(x, a) / tgamma(a + 1);
} else {
result = Math.pow(x, a) / a;
}
result *= 1 - a * x / (a + 1);
break;
case 7:
default:
// x is large,
// Compute Q:
result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
result /= x;
if (result != 0) {
result *= incompleteTgammaLargeX(a, x, pol);
}
break;
}
if (normalised && (result > 1)) {
result = 1;
}
if (invert) {
final double gam = normalised ? 1 : tgamma(a);
result = gam - result;
}
return result;
}