You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/12/13 00:14:31 UTC
[commons-numbers] branch master updated: Move the direct computation before checks for base near 1
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
The following commit(s) were added to refs/heads/master by this push:
new 1084917 Move the direct computation before checks for base near 1
1084917 is described below
commit 10849175095e3569e2f096f3f967ce699e39be0a
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Mon Dec 13 00:10:54 2021 +0000
Move the direct computation before checks for base near 1
When the direct computation can be used then there is no requirement for
the Lanczos code that uses pow(z / (z+g-0.5), a). Thus there is no need
to avoid a base near 1 in the power function. This check only applies
after the direct computation has been eliminated.
---
.../apache/commons/numbers/gamma/BoostGamma.java | 43 +++++++++++-----------
1 file changed, 22 insertions(+), 21 deletions(-)
diff --git a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java
index 76feed6..0a5b9a1 100644
--- a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java
+++ b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostGamma.java
@@ -1794,6 +1794,27 @@ final class BoostGamma {
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 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
@@ -1803,6 +1824,7 @@ final class BoostGamma {
double prefix;
final double d = ((z - a) - Lanczos.gmh()) / agh;
+
if ((Math.abs(d * d * a) <= 100) && (a > 150)) {
// 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:
@@ -1811,27 +1833,6 @@ final class BoostGamma {
prefix = a * SpecialMath.log1pmx(d) + z * -Lanczos.gmh() / agh;
prefix = Math.exp(prefix);
} else {
- // 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 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);
- }
- }
-
//
// general case.
// direct computation is most accurate, but use various fallbacks