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
commit ab1cc9e7ed7b471bd4d61578ff0d08b1a5c67520
Author: aherbert <aherb...@apache.org>
AuthorDate: Tue Aug 29 12:37:36 2023 +0100
Replace extended precision methods in gamma erf functions with DD
The computation of the round-off error from squaring a number can be
performed using DD.ofSquare(double).
---
.../org/apache/commons/numbers/gamma/BoostErf.java | 102 ++-------------------
1 file changed, 6 insertions(+), 96 deletions(-)
diff --git
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
index b95e3526..476c6562 100644
---
a/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
+++
b/commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/BoostErf.java
@@ -22,6 +22,8 @@
package org.apache.commons.numbers.gamma;
+import org.apache.commons.numbers.core.DD;
+
/**
* Implementation of the <a href="http://mathworld.wolfram.com/Erf.html";>error
function</a> and
* its inverse.
@@ -39,13 +41,6 @@ package org.apache.commons.numbers.gamma;
* Boost C++ Error functions</a>
*/
final class BoostErf {
- /**
- * The multiplier used to split the double value into high and low parts.
From
- * Dekker (1971): "The constant should be chosen equal to 2^(p - p/2) + 1,
- * where p is the number of binary digits in the mantissa". Here p is 53
- * and the multiplier is {@code 2^27 + 1}.
- */
- private static final double MULTIPLIER = 1.0 + 0x1.0p27;
/** 1 / sqrt(pi). Used for the scaled complementary error function erfcx.
*/
private static final double ONE_OVER_ROOT_PI =
0.5641895835477562869480794515607725858;
/** Threshold for the scaled complementary error function erfcx
@@ -734,9 +729,8 @@ final class BoostErf {
// round-off b is negative or zero:
// exp(a) * exp1m(b) + exp(a)
// inf * 0 + inf or inf * -b + inf
- final double a = x * x;
- final double b = squareLowUnscaled(x, a);
- return expxx(a, b);
+ final DD x2 = DD.ofSquare(x);
+ return expxx(x2.hi(), x2.lo());
}
/**
@@ -754,9 +748,8 @@ final class BoostErf {
* @return exp(-x*x)
*/
static double expmxx(double x) {
- final double a = x * x;
- final double b = squareLowUnscaled(x, a);
- return expxx(-a, -b);
+ final DD x2 = DD.ofSquare(x);
+ return expxx(-x2.hi(), -x2.lo());
}
/**
@@ -792,87 +785,4 @@ final class BoostErf {
// b ~ expm1(b)
return ea * b + ea;
}
-
- // Extended precision multiplication specialised for the square adapted
from:
- // org.apache.commons.numbers.core.ExtendedPrecision
-
- /**
- * Compute the low part of the double length number {@code (z,zz)} for the
exact
- * square of {@code x} using Dekker's mult12 algorithm. The standard
precision product
- * {@code x*x} must be provided. The number {@code x} is split into high
and low parts
- * using Dekker's algorithm.
- *
- * <p>Warning: This method does not perform scaling in Dekker's split and
large
- * finite numbers can create NaN results.
- *
- * @param x Number to square
- * @param xx Standard precision product {@code x*x}
- * @return the low part of the square double length number
- */
- private static double squareLowUnscaled(double x, double xx) {
- // Split the numbers using Dekker's algorithm without scaling
- final double hx = highPartUnscaled(x);
- final double lx = x - hx;
-
- return squareLow(hx, lx, xx);
- }
-
- /**
- * Implement Dekker's method to split a value into two parts. Multiplying
by (2^s + 1) creates
- * a big value from which to derive the two split parts.
- * <pre>
- * c = (2^s + 1) * a
- * a_big = c - a
- * a_hi = c - a_big
- * a_lo = a - a_hi
- * a = a_hi + a_lo
- * </pre>
- *
- * <p>The multiplicand allows a p-bit value to be split into
- * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value
{@code a_lo}.
- * Combined they have (p-1) bits of significand but the sign bit of {@code
a_lo}
- * contains a bit of information. The constant is chosen so that s is
ceil(p/2) where
- * the precision p for a double is 53-bits (1-bit of the mantissa is
assumed to be
- * 1 for a non sub-normal number) and s is 27.
- *
- * <p>This conversion does not use scaling and the result of overflow is
NaN. Overflow
- * may occur when the exponent of the input value is above 996.
- *
- * <p>Splitting a NaN or infinite value will return NaN.
- *
- * @param value Value.
- * @return the high part of the value.
- * @see Math#getExponent(double)
- */
- private static double highPartUnscaled(double value) {
- final double c = MULTIPLIER * value;
- return c - (c - value);
- }
-
- /**
- * Compute the low part of the double length number {@code (z,zz)} for the
exact
- * square of {@code x} using Dekker's mult12 algorithm. The standard
- * precision product {@code x*x} must be provided. The number {@code x}
- * should already be split into low and high parts.
- *
- * <p>Note: This uses the high part of the result {@code (z,zz)} as {@code
x * x} and not
- * {@code hx * hx + hx * lx + lx * hx} as specified in Dekker's original
paper.
- * See Shewchuk (1997) for working examples.
- *
- * @param hx High part of factor.
- * @param lx Low part of factor.
- * @param xx Square of the factor.
- * @return <code>lx * ly - (((xy - hx * hy) - lx * hy) - hx * ly)</code>
- * @see <a
href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps";>
- * Shewchuk (1997) Theorum 18</a>
- */
- private static double squareLow(double hx, double lx, double xx) {
- // Compute the multiply low part:
- // err1 = xy - hx * hy
- // err2 = err1 - lx * hy
- // err3 = err2 - hx * ly
- // low = lx * ly - err3
- return lx * lx - ((xx - hx * hx) - 2 * lx * hx);
- }
}
-
