in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostBeta.java [255:397]
private static double ibetaPowerTerms(double a, double b, double x,
double y, boolean normalised, double prefix) {
if (!normalised) {
// can we do better here?
return Math.pow(x, a) * Math.pow(y, b);
}
double result;
final double c = a + b;
// combine power terms with Lanczos approximation:
final double agh = a + BoostGamma.Lanczos.GMH;
final double bgh = b + BoostGamma.Lanczos.GMH;
final double cgh = c + BoostGamma.Lanczos.GMH;
result = BoostGamma.Lanczos.lanczosSumExpGScaled(c) /
(BoostGamma.Lanczos.lanczosSumExpGScaled(a) * BoostGamma.Lanczos.lanczosSumExpGScaled(b));
result *= prefix;
// combine with the leftover terms from the Lanczos approximation:
result *= Math.sqrt(bgh / Math.E);
result *= Math.sqrt(agh / cgh);
// l1 and l2 are the base of the exponents minus one:
double l1 = (x * b - y * agh) / agh;
double l2 = (y * a - x * bgh) / bgh;
if (Math.min(Math.abs(l1), Math.abs(l2)) < 0.2) {
// when the base of the exponent is very near 1 we get really
// gross errors unless extra care is taken:
if (l1 * l2 > 0 || Math.min(a, b) < 1) {
//
// This first branch handles the simple cases where either:
//
// * The two power terms both go in the same direction
// (towards zero or towards infinity). In this case if either
// term overflows or underflows, then the product of the two must
// do so also.
// * Alternatively if one exponent is less than one, then we
// can't productively use it to eliminate overflow or underflow
// from the other term. Problems with spurious overflow/underflow
// can't be ruled out in this case, but it is *very* unlikely
// since one of the power terms will evaluate to a number close to 1.
//
if (Math.abs(l1) < 0.1) {
result *= Math.exp(a * Math.log1p(l1));
} else {
result *= Math.pow((x * cgh) / agh, a);
}
if (Math.abs(l2) < 0.1) {
result *= Math.exp(b * Math.log1p(l2));
} else {
result *= Math.pow((y * cgh) / bgh, b);
}
} else if (Math.max(Math.abs(l1), Math.abs(l2)) < 0.5) {
//
// Both exponents are near one and both the exponents are
// greater than one and further these two
// power terms tend in opposite directions (one towards zero,
// the other towards infinity), so we have to combine the terms
// to avoid any risk of overflow or underflow.
//
// We do this by moving one power term inside the other, we have:
//
// (1 + l1)^a * (1 + l2)^b
// = ((1 + l1)*(1 + l2)^(b/a))^a
// = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1
// = exp((b/a) * log(1 + l2)) - 1
//
// The tricky bit is deciding which term to move inside :-)
// By preference we move the larger term inside, so that the
// size of the largest exponent is reduced. However, that can
// only be done as long as l3 (see above) is also small.
//
final boolean smallA = a < b;
final double ratio = b / a;
if ((smallA && ratio * l2 < 0.1) || (!smallA && l1 / ratio > 0.1)) {
double l3 = Math.expm1(ratio * Math.log1p(l2));
l3 = l1 + l3 + l3 * l1;
l3 = a * Math.log1p(l3);
result *= Math.exp(l3);
} else {
double l3 = Math.expm1(Math.log1p(l1) / ratio);
l3 = l2 + l3 + l3 * l2;
l3 = b * Math.log1p(l3);
result *= Math.exp(l3);
}
} else if (Math.abs(l1) < Math.abs(l2)) {
// First base near 1 only:
double l = a * Math.log1p(l1) + b * Math.log((y * cgh) / bgh);
if (l <= LOG_MIN_VALUE || l >= LOG_MAX_VALUE) {
l += Math.log(result);
// Update: Allow overflow to infinity
result = Math.exp(l);
} else {
result *= Math.exp(l);
}
} else {
// Second base near 1 only:
double l = b * Math.log1p(l2) + a * Math.log((x * cgh) / agh);
if (l <= LOG_MIN_VALUE || l >= LOG_MAX_VALUE) {
l += Math.log(result);
// Update: Allow overflow to infinity
result = Math.exp(l);
} else {
result *= Math.exp(l);
}
}
} else {
// general case:
final double b1 = (x * cgh) / agh;
final double b2 = (y * cgh) / bgh;
l1 = a * Math.log(b1);
l2 = b * Math.log(b2);
if (l1 >= LOG_MAX_VALUE || l1 <= LOG_MIN_VALUE || l2 >= LOG_MAX_VALUE || l2 <= LOG_MIN_VALUE) {
// Oops, under/overflow, sidestep if we can:
if (a < b) {
final double p1 = Math.pow(b2, b / a);
final double l3 = a * (Math.log(b1) + Math.log(p1));
if (l3 < LOG_MAX_VALUE && l3 > LOG_MIN_VALUE) {
result *= Math.pow(p1 * b1, a);
} else {
l2 += l1 + Math.log(result);
// Update: Allow overflow to infinity
result = Math.exp(l2);
}
} else {
final double p1 = Math.pow(b1, a / b);
final double l3 = (Math.log(p1) + Math.log(b2)) * b;
if (l3 < LOG_MAX_VALUE && l3 > LOG_MIN_VALUE) {
result *= Math.pow(p1 * b2, b);
} else {
l2 += l1 + Math.log(result);
// Update: Allow overflow to infinity
result = Math.exp(l2);
}
}
} else {
// finally the normal case:
result *= Math.pow(b1, a) * Math.pow(b2, b);
}
}
return result;
}