This is an automated email from the ASF dual-hosted git repository.
asf-gitbox-commits 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 bcbc8f12 Add reference error test for DD.pow
bcbc8f12 is described below
commit bcbc8f12713dc7867d36b78c6f56a1b752241c44
Author: Alex Herbert <[email protected]>
AuthorDate: Mon Oct 19 23:52:39 2026 +0100
Add reference error test for DD.pow
Add alternative pow implementation to DDExt test class for comparison.
The current implementation has lower RMS and max error.
---
.../org/apache/commons/numbers/core/DDExt.java | 166 +++++++++++++++++++++
.../org/apache/commons/numbers/core/DDTest.java | 88 ++++++++++-
.../org/apache/commons/numbers/core/TestUtils.java | 4 +-
3 files changed, 254 insertions(+), 4 deletions(-)
diff --git
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDExt.java
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDExt.java
index 133237c1..c3bbddd4 100644
---
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDExt.java
+++
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDExt.java
@@ -732,4 +732,170 @@ public final class DDExt {
exp[0] = b * n + qb + rb + expi[0];
return f;
}
+
+ /**
+ * Compute number {@code x} raised to the power {@code n}.
+ *
+ * <p>Note: This is an alternative implementation to {@link DD#pow(int,
long[])}.
+ * Accuracy is lower (see {@link DDTest#computePowError()}.
+ *
+ * <p>The value is returned as fractional {@code f} and integral
+ * {@code 2^exp} components.
+ * <pre>
+ * (x+xx)^n = (f+ff) * 2^exp
+ * </pre>
+ *
+ * <p>The combined fractional part (f, ff) is in the range {@code [0.5,
1)}.
+ *
+ * <p>Special cases:
+ * <ul>
+ * <li>If {@code (x, xx)} is zero the high part of the fractional part is
+ * computed using {@link Math#pow(double, double) Math.pow(x, n)} and
the exponent is 0.</li>
+ * <li>If {@code n = 0} the fractional part is 0.5 and the exponent is
1.</li>
+ * <li>If {@code (x, xx)} is an exact power of 2 the fractional part is
0.5 and the exponent
+ * is the power of 2 minus 1.</li>
+ * <li>If the result high-part is an exact power of 2 and the low-part
has an opposite
+ * signed non-zero magnitude then the fraction high-part {@code f}
will be {@code +/-1} such that
+ * the double-double number is in the range {@code [0.5, 1)}.</li>
+ * <li>If the argument is not finite then a fractional representation is
not possible.
+ * In this case the fraction and the scale factor is undefined.</li>
+ * </ul>
+ *
+ * <p>The computed result is approximately {@code 16 * (n - 1) * eps} of
the exact result
+ * where eps is 2<sup>-106</sup>.
+ *
+ * @param x Number.
+ * @param n Power.
+ * @param exp Result power of two scale factor (integral exponent).
+ * @return Fraction part.
+ * @see DD#frexp(int[])
+ */
+ public static DD pow(DD x, int n, long[] exp) {
+ // Edge cases.
+ if (n == 0) {
+ exp[0] = 1;
+ return DD.of(0.5, 0);
+ }
+ // IEEE result for non-finite or zero
+ if (!Double.isFinite(x.hi()) || x.hi() == 0) {
+ exp[0] = 0;
+ return DD.of(Math.pow(x.hi(), n), 0);
+ }
+ // Here the number is non-zero finite
+ final int[] ie = {0};
+ DD f = x.frexp(ie);
+ final long b = ie[0];
+ // Handle exact powers of 2
+ if (Math.abs(f.hi()) == 0.5 && f.lo() == 0) {
+ // (f * 2^b)^n = (2f)^n * 2^(b-1)^n
+ // Use Math.pow to create the sign.
+ // Note the result must be scaled to the fractional representation
+ // by multiplication by 0.5 and addition of 1 to the exponent.
+ final double y0 = 0.5 * Math.pow(2 * f.hi(), n);
+ // Propagate sign change (y0*f.hi()) to the original zero
(this.lo())
+ final double y1 = Math.copySign(0.0, y0 * f.hi() * x.lo());
+ exp[0] = 1 + (b - 1) * n;
+ return DD.of(y0, y1);
+ }
+ if (n < 0) {
+ f = computePowScaled(b, f.hi(), f.lo(), -n, exp);
+ // Result is a non-zero fraction part so inversion is safe
+ f = DD.of(f.hi(), f.lo()).reciprocal();
+ // Rescale to [0.5, 1.0)
+ f = f.frexp(ie);
+ exp[0] = ie[0] - exp[0];
+ return f;
+ }
+ return computePowScaled(b, f.hi(), f.lo(), n, exp);
+ }
+
+ /**
+ * Compute the number {@code x} (non-zero finite) raised to the power
{@code n}.
+ *
+ * <p>The input power is treated as an unsigned integer. Thus the negative
value
+ * {@link Integer#MIN_VALUE} is 2^31.
+ *
+ * @param b Integral component 2^exp of x.
+ * @param x Fractional high part of x.
+ * @param xx Fractional low part of x.
+ * @param n Power (in [2, 2^31]).
+ * @param exp Result power of two scale factor (integral exponent).
+ * @return Fraction part.
+ */
+ private static DD computePowScaled(long b, double x, double xx, int n,
long[] exp) {
+ // Compute the power by multiplication (keeping track of the scale):
+ // 13 = 1101
+ // x^13 = x^8 * x^4 * x^1
+ //
+ // 21 = 10101
+ // x^21 = x^16 * x^4 * x^1
+ //
+ // 1. Initialise result r as 1
+ // 2. Initialise exponent nn as n
+ // 3. Initialise square x2 as x
+ // 4. while remaining bits (0 or 1):
+ // - if the lowest bit is 1 then multiply r by x2
+ // - consume bit and if none left stop
+ // - square x2
+
+ long re = 0;
+ DD r = DD.ONE;
+
+ // Initialise the square as x^1. Represented as 2^fe * f.
+ // Scale the input in [0.5, 1) to be above 1 so multiplication
increases magnitude.
+ long fe = b - 1;
+ double f0 = x * 2;
+ double f1 = xx * 2;
+
+ double u;
+ double v;
+ double w;
+
+ int nn = n;
+ while (true) {
+ if ((nn & 0x1) != 0) {
+ // r = r * x2
+ r = r.multiply(DD.of(f0, f1));
+ re += fe;
+ if (Math.abs(r.hi()) > SAFE_MULTIPLY) {
+ // Scale back to the [1, 2) range. As safe multiply is
2^500
+ // the exponent should be < 1001 so the twoPow scaling
factor is supported.
+ final int e = Math.getExponent(r.hi());
+ final double s = DD.twoPow(-e);
+ re += e;
+ r = DD.of(r.hi() * s, r.lo() * s);
+ }
+
+ }
+ nn >>>= 1;
+ if (nn == 0) {
+ break;
+ }
+ // Square x2
+ // Inline multiply(f0, f1, f0, f1), adapted for squaring
+ fe <<= 1;
+ u = f0 * f0;
+ v = DD.twoSquareLow(f0, u);
+ // Inline fastTwoSum(hi, lo + (2 * f0 * f1), f)
+ w = v + (2 * f0 * f1);
+ f0 = u + w;
+ f1 = DD.fastTwoSumLow(u, w, f0);
+ // Rescale
+ if (Math.abs(f0) > SAFE_MULTIPLY) {
+ // Scale back to the [1, 2) range. As safe multiply is 2^500
+ // the exponent should be < 1001 so the twoPow scaling factor
is supported.
+ final int e = Math.getExponent(f0);
+ final double s = DD.twoPow(-e);
+ fe += e;
+ f0 *= s;
+ f1 *= s;
+ }
+ }
+
+ // Normalise the return
+ final int[] e = {0};
+ r = r.frexp(e);
+ exp[0] = re + e[0];
+ return r;
+ }
}
diff --git
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDTest.java
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDTest.java
index 06b047e2..467b93ac 100644
---
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDTest.java
+++
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DDTest.java
@@ -3033,6 +3033,29 @@ class DDTest {
// a -> a.abs().sqrt(MC_DIVIDE), a -> a.abs().sqrt(), a ->
DDExt.sqrt(a.abs()));
}
+ /**
+ * Compute the max and mean error of scaled power operations.
+ * Note: This is not a test and is included for reference.
+ */
+ @Test
+ @org.junit.jupiter.api.Disabled("This fixture creates reference data")
+ void computePowError() {
+ final int n = 10_000;
+ final int samples = 1_000_000;
+ final UniformRandomProvider rng = createRNG();
+ final DD[] dd = IntStream.range(0, n).mapToObj(i ->
+ signedNormalDoubleDouble(rng)
+ ).toArray(DD[]::new);
+ final BigDecimal[] bd =
Arrays.stream(dd).map(DD::bigDecimalValue).toArray(BigDecimal[]::new);
+ // BigDecimal power supported up to 2^(2^29)
+ for (int b = 2; b < 29; b++) {
+ computePowError(rng, dd, bd, Math.min(n, samples), b);
+ }
+ for (int b = 2; b < 29; b++) {
+ computePowError(rng, dd, bd, Math.min(n, samples), -b);
+ }
+ }
+
/**
* Compute the max and mean error of the operators. This method supports
two operators
* so they can be compared to a reference implementation.
@@ -3042,7 +3065,7 @@ class DDTest {
* @param bd DD data converted to BigDecimal
* @param samples Number of samples
* @param operatorName Operator name
- * @param expected Operator to compute the expected result
+ * @param reference Operator to compute the expected result
* @param op1 First operator
* @param op2 Second operator
*/
@@ -3083,7 +3106,7 @@ class DDTest {
* @param bd DD data converted to BigDecimal
* @param samples Number of samples
* @param operatorName Operator name
- * @param expected Operator to compute the expected result
+ * @param reference Operator to compute the expected result
* @param op1 First operator
* @param op2 Second operator
*/
@@ -3116,6 +3139,67 @@ class DDTest {
1e-9 * time);
}
+ /**
+ * Compute the max and mean error of the pow operator.
+ * The javadoc states 16 * (n - 1) eps. This is an upper error bound.
+ * This reports RMS error ~ n/4; Max error ~ 2n.
+ *
+ * <p>The alternative implementation is {@link DDExt#pow(DD, int,
long[])}. This computes
+ * the result starting at 1 and multiplying by increasing powers of 2 if
required. The error is
+ * worse than the DD power implementation which starts at x^1 and
multiplies by x if required,
+ * squaring the current result each iteration. Due to pre-computation of
factors to multiply
+ * by x, the current DD implementation is also faster.
+ *
+ * <p>This method can be used to show DDMath.pow is accurate to ~ 0.5 eps.
+ *
+ * @param rng Source of randomness
+ * @param dd DD data
+ * @param bd DD data converted to BigDecimal
+ * @param samples Number of samples
+ * @param b Maximum bit in the power
+ */
+ private static void computePowError(UniformRandomProvider rng,
+ DD[] dd, BigDecimal[] bd, int samples, int b) {
+ final int absB = Math.abs(b);
+ Assertions.assertTrue(-absB >= -30, "b is not an exponent; exponent is
2^b");
+ final ErrorStatistics e1 = new ErrorStatistics();
+ final ErrorStatistics e2 = new ErrorStatistics();
+ final long start = System.nanoTime();
+ final int[] count = {0};
+ final int base = 1 << absB;
+ final IntSupplier next = samples == dd.length ? () -> count[0]++ : ()
-> rng.nextInt(dd.length);
+ final long[] exp = {0};
+ final MathContext mc = MathContext.DECIMAL128;
+ final BigDecimal two = BigDecimal.valueOf(2);
+ final BigDecimal half = new BigDecimal(0.5);
+ for (int n = 0; n < samples; n++) {
+ final int i = next.getAsInt();
+ int y = base + rng.nextInt(base);
+ if (b < 0) {
+ y = -y;
+ }
+ BigDecimal expected = bd[i].pow(y, mc);
+ final DD actual1 = dd[i].pow(y, exp);
+ e1.add(TestUtils.assertScaledEquals(expected, actual1, exp[0],
-1.0, () -> "DD.pow"));
+ final DD actual2 = DDExt.pow(dd[i], y, exp);
+ // Max ULP error ~ 0.55:
+ // final DD actual2 = DDMath.pow(dd[i], y, exp);
+ e2.add(TestUtils.assertScaledEquals(expected, actual2, exp[0],
-1.0, () -> "DDExt.pow"));
+ }
+ final long time = System.nanoTime() - start;
+ final double rms1 = e1.getRMS() / EPS;
+ final double max1 = e1.getMaxAbs() / EPS;
+ final double rms2 = e2.getRMS() / EPS;
+ final double max2 = e2.getMaxAbs() / EPS;
+ // Report errors scaled by the maximum power
+ final double scale = Math.scalb(1.0, -absB - 1);
+ TestUtils.printf("pow 2^%-3d %12.6g (%4.2f) %12.6g (%4.2f) %12.6g
%12.6g (%4.2f) %12.6g (%4.2f) %12.6g (%.3fs)%n",
+ b,
+ rms1, scale * rms1, max1, scale * max1, e1.getMean() / EPS,
+ rms2, scale * rms2, max2, scale * max2, e2.getMean() / EPS,
+ 1e-9 * time);
+ }
+
/**
* Class to compute the error statistics.
*
diff --git
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/TestUtils.java
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/TestUtils.java
index f936a9c8..8615308e 100644
---
a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/TestUtils.java
+++
b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/TestUtils.java
@@ -218,9 +218,9 @@ public final class TestUtils {
// e * 2^exp = expected
// e = expected * 2^-exp
if (exp < 0) {
- e = expected.multiply(TWO.pow((int) -exp), MC_SCALING);
+ e = expected.multiply(TWO.pow((int) -exp, MC_SCALING), MC_SCALING);
} else {
- e = expected.multiply(HALF.pow((int) exp), MC_SCALING);
+ e = expected.multiply(HALF.pow((int) exp, MC_SCALING), MC_SCALING);
}
return TestUtils.assertEquals(e, actual, eps, () -> prefix(msg) +
"scale=2^" + exp);
}
