in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java [1716:1822]
static double regularisedGammaPrefix(double a, double z) {
if (z >= Double.MAX_VALUE) {
return 0;
}
// Update this condition from: a < 1
if (a <= 1) {
//
// We have to treat a < 1 as a special case because our Lanczos
// approximations are optimised against the factorials with a > 1,
// and for high precision types especially (128-bit reals for example)
// very small values of a can give rather erroneous results for gamma
// unless we do this:
//
// Boost todo: is this still required? Lanczos approx should be better now?
//
// Update this condition from: z <= LOG_MIN_VALUE
// Require exp(-z) to not underflow:
// -z > log(min_value)
if (-z <= LOG_MIN_VALUE) {
// Oh dear, have to use logs, should be free of cancellation errors though:
return Math.exp(a * Math.log(z) - z - lgamma(a));
}
// direct calculation, no danger of overflow as gamma(a) < 1/a
// for small a.
return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
}
// Update to the Boost code.
// Use some of the logic from fullIgammaPrefix(a, z) to use the direct
// computation if it is valid. Assuming pow and exp are accurate to 1 ULP it
// puts most of the error in evaluation of tgamma(a). This is accurate
// enough that this reduces max error on the current test data.
//
// Overflow cases fall-through to the Lanczos approximation that incorporates
// the pow and exp terms used in the tgamma(a) computation with the terms
// z^a and e^-z into a single evaluation of pow and exp. See equation 15:
// https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
if (a <= MAX_GAMMA_Z) {
final double alz1 = a * Math.log(z);
if (z >= 1) {
if ((alz1 < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
}
} else if (alz1 > LOG_MIN_VALUE) {
return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
}
}
//
// For smallish a and x combining the power terms with the Lanczos approximation
// gives the greatest accuracy
//
final double agh = a + Lanczos.GMH;
double prefix;
final double factor = Math.sqrt(agh / Math.E) / Lanczos.lanczosSumExpGScaled(a);
// Update to the Boost code.
// Lower threshold for large a from 150 to 128 and compute d on demand.
// See NUMBERS-179.
if (a > 128) {
final double d = ((z - a) - Lanczos.GMH) / agh;
if (Math.abs(d * d * a) <= 100) {
// special case for large a and a ~ z.
// When a and x are large, we end up with a very large exponent with a base near one:
// this will not be computed accurately via the pow function, and taking logs simply
// leads to cancellation errors.
prefix = a * SpecialMath.log1pmx(d) + z * -Lanczos.GMH / agh;
prefix = Math.exp(prefix);
return prefix * factor;
}
}
//
// general case.
// direct computation is most accurate, but use various fallbacks
// for different parts of the problem domain:
//
final double alz = a * Math.log(z / agh);
final double amz = a - z;
if ((Math.min(alz, amz) <= LOG_MIN_VALUE) || (Math.max(alz, amz) >= LOG_MAX_VALUE)) {
final double amza = amz / a;
if ((Math.min(alz, amz) / 2 > LOG_MIN_VALUE) && (Math.max(alz, amz) / 2 < LOG_MAX_VALUE)) {
// compute square root of the result and then square it:
final double sq = Math.pow(z / agh, a / 2) * Math.exp(amz / 2);
prefix = sq * sq;
} else if ((Math.min(alz, amz) / 4 > LOG_MIN_VALUE) &&
(Math.max(alz, amz) / 4 < LOG_MAX_VALUE) && (z > a)) {
// compute the 4th root of the result then square it twice:
final double sq = Math.pow(z / agh, a / 4) * Math.exp(amz / 4);
prefix = sq * sq;
prefix *= prefix;
} else if ((amza > LOG_MIN_VALUE) && (amza < LOG_MAX_VALUE)) {
prefix = Math.pow((z * Math.exp(amza)) / agh, a);
} else {
prefix = Math.exp(alz + amz);
}
} else {
prefix = Math.pow(z / agh, a) * Math.exp(amz);
}
prefix *= factor;
return prefix;
}