Update threshold in regularisedGammaPrefix
Modifies the threshold from a > 150 (Boost version) to a > 128.
@@ -1822,16 +1822,25 @@ final class BoostGamma {
final double agh = a + Lanczos.gmh();
double prefix;
- final double d = ((z - a) - Lanczos.gmh()) / agh;
- if ((Math.abs(d * d * a) <= 100) && (a > 100)) {
+ 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);
- } else {
+ return prefix * factor;
+ }
+ }
+
//
// general case.
// direct computation is most accurate, but use various fallbacks
@@ -1860,8 +1869,7 @@ final class BoostGamma {
} else {
prefix = Math.pow(z / agh, a) * Math.exp(amz);
}
- }
- prefix *= Math.sqrt(agh / Math.E) / Lanczos.lanczosSumExpGScaled(a);
+ prefix *= factor;
return prefix;
}
