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 644347bd Allow use of different multiplication in DDExt.pow
644347bd is described below
commit 644347bdad17698d9ab724602c06134d9454340d
Author: Alex Herbert <[email protected]>
AuthorDate: Tue Oct 20 11:17:16 2026 +0100
Allow use of different multiplication in DDExt.pow
---
.../org/apache/commons/numbers/core/DDExt.java | 78 +++++++++++++---------
.../org/apache/commons/numbers/core/DDTest.java | 49 ++++++++++----
2 files changed, 80 insertions(+), 47 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 c3bbddd4..8c7205d2 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
@@ -16,6 +16,9 @@
*/
package org.apache.commons.numbers.core;
+import java.util.function.BiFunction;
+import java.util.function.Function;
+
/**
* Computes double-double floating-point operations.
*
@@ -737,7 +740,19 @@ public final class DDExt {
* 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()}.
+ * It initialises the result as 1 and selectively multiplies by an
increasing magnitude term
+ * which is squared each iteration.
+ *
+ * <p>The method used for multiplication and square can be specified.
+ * <ol>
+ * <li>{@link DD#multiply(DD)}
+ * <li>{@link DDExt#multiply(DD, DD)}
+ * </ol>
+ *
+ * <p>When using DD multiplication the accuracy is lower than the DD pow
implementation
+ * (see {@link DDTest#computePowError()}.
+ * When using accurate multiplication from this class accuracy improves
approximately
+ * 4-fold for a performance penalty of 4.5-fold. The method is not ulp
accurate.
*
* <p>The value is returned as fractional {@code f} and integral
* {@code 2^exp} components.
@@ -761,16 +776,14 @@ public final class DDExt {
* 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).
+ * @param mode Multiplication mode.
* @return Fraction part.
* @see DD#frexp(int[])
*/
- public static DD pow(DD x, int n, long[] exp) {
+ public static DD pow(DD x, int n, long[] exp, int mode) {
// Edge cases.
if (n == 0) {
exp[0] = 1;
@@ -798,7 +811,7 @@ public final class DDExt {
return DD.of(y0, y1);
}
if (n < 0) {
- f = computePowScaled(b, f.hi(), f.lo(), -n, exp);
+ f = computePowScaled(b, f.hi(), f.lo(), -n, exp, mode);
// 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)
@@ -806,7 +819,7 @@ public final class DDExt {
exp[0] = ie[0] - exp[0];
return f;
}
- return computePowScaled(b, f.hi(), f.lo(), n, exp);
+ return computePowScaled(b, f.hi(), f.lo(), n, exp, mode);
}
/**
@@ -820,9 +833,10 @@ public final class DDExt {
* @param xx Fractional low part of x.
* @param n Power (in [2, 2^31]).
* @param exp Result power of two scale factor (integral exponent).
+ * @param mode Multiplication mode.
* @return Fraction part.
*/
- private static DD computePowScaled(long b, double x, double xx, int n,
long[] exp) {
+ private static DD computePowScaled(long b, double x, double xx, int n,
long[] exp, int mode) {
// Compute the power by multiplication (keeping track of the scale):
// 13 = 1101
// x^13 = x^8 * x^4 * x^1
@@ -838,25 +852,31 @@ public final class DDExt {
// - consume bit and if none left stop
// - square x2
+ // Note: This implementation is not optimised. It is used to test use
+ // of different multiplication routines for the power result.
+ final BiFunction<DD, DD, DD> multiply;
+ final Function<DD, DD> square;
+ if (mode == 1) {
+ multiply = DD::multiply;
+ square = DD::square;
+ } else {
+ multiply = DDExt::multiply;
+ square = DDExt::square;
+ }
+
long re = 0;
DD r = DD.ONE;
- // Initialise the square as x^1. Represented as 2^fe * f.
+ // Initialise the square as x.
// 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;
+ long x2e = b - 1;
+ DD x2 = DD.of(x * 2, xx * 2);
int nn = n;
while (true) {
if ((nn & 0x1) != 0) {
- // r = r * x2
- r = r.multiply(DD.of(f0, f1));
- re += fe;
+ r = multiply.apply(r, x2);
+ re += x2e;
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.
@@ -871,24 +891,16 @@ public final class DDExt {
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);
+ x2e <<= 1;
+ x2 = square.apply(x2);
// Rescale
- if (Math.abs(f0) > SAFE_MULTIPLY) {
+ if (Math.abs(x2.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(f0);
+ final int e = Math.getExponent(x2.hi());
final double s = DD.twoPow(-e);
- fe += e;
- f0 *= s;
- f1 *= s;
+ x2e += e;
+ x2 = DD.of(x2.hi() * s, x2.lo() * s);
}
}
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 467b93ac..ae18f545 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
@@ -22,6 +22,7 @@ import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.Arrays;
+import java.util.function.BiFunction;
import java.util.function.BinaryOperator;
import java.util.function.IntSupplier;
import java.util.function.Supplier;
@@ -3047,12 +3048,14 @@ class DDTest {
signedNormalDoubleDouble(rng)
).toArray(DD[]::new);
final BigDecimal[] bd =
Arrays.stream(dd).map(DD::bigDecimalValue).toArray(BigDecimal[]::new);
+ // Use to choose the alternative implementation
+ final int alternative = 1;
// BigDecimal power supported up to 2^(2^29)
for (int b = 2; b < 29; b++) {
- computePowError(rng, dd, bd, Math.min(n, samples), b);
+ computePowError(rng, dd, bd, Math.min(n, samples), b, alternative);
}
for (int b = 2; b < 29; b++) {
- computePowError(rng, dd, bd, Math.min(n, samples), -b);
+ computePowError(rng, dd, bd, Math.min(n, samples), -b,
alternative);
}
}
@@ -3141,14 +3144,22 @@ class DDTest {
/**
* Compute the max and mean error of the pow operator.
- * The javadoc states 16 * (n - 1) eps. This is an upper error bound.
+ *
+ * <p>The DD.pow 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>The alternative implementation is a choice of:
+ * <ol>
+ * <li>{@link DDExt#pow(DD, int, long[])}
+ * <li>{@link DDExt#accuratePow(DD, int, long[])}
+ * <li>{@link DDMath#pow(DD, int, long[])}
+ * </ol>
+ *
+ * <p>Note that [1] and [2] compute the result starting at 1 and
multiplying by
+ * increasing powers of 2 if required. In the case of [1] 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.
*
@@ -3157,9 +3168,10 @@ class DDTest {
* @param bd DD data converted to BigDecimal
* @param samples Number of samples
* @param b Maximum bit in the power
+ * @param alternative Alternative power implementation.
*/
private static void computePowError(UniformRandomProvider rng,
- DD[] dd, BigDecimal[] bd, int samples, int b) {
+ DD[] dd, BigDecimal[] bd, int samples, int b, int alternative) {
final int absB = Math.abs(b);
Assertions.assertTrue(-absB >= -30, "b is not an exponent; exponent is
2^b");
final ErrorStatistics e1 = new ErrorStatistics();
@@ -3172,6 +3184,17 @@ class DDTest {
final MathContext mc = MathContext.DECIMAL128;
final BigDecimal two = BigDecimal.valueOf(2);
final BigDecimal half = new BigDecimal(0.5);
+ // The alternative implementation uses the exp variable for output
+ final BiFunction<DD, Integer, DD> pow2;
+ final String name;
+ if (alternative == 3) {
+ pow2 = (x, n) -> DDMath.pow(x, n, exp);
+ name = "DDMath.pow";
+ } else {
+ // Supports alternative 1 or 2
+ pow2 = (x, n) -> DDExt.pow(x, n, exp, alternative);
+ name = alternative == 1 ? "DDExt.pow" : "DDExt.accPow";
+ }
for (int n = 0; n < samples; n++) {
final int i = next.getAsInt();
int y = base + rng.nextInt(base);
@@ -3181,9 +3204,7 @@ class DDTest {
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);
+ final DD actual2 = pow2.apply(dd[i], y);
e2.add(TestUtils.assertScaledEquals(expected, actual2, exp[0],
-1.0, () -> "DDExt.pow"));
}
final long time = System.nanoTime() - start;
@@ -3193,8 +3214,8 @@ class DDTest {
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,
+ TestUtils.printf("%-13s 2^%-3d %12.6g (%4.2f) %12.6g (%4.2f) %12.6g
%12.6g (%4.2f) %12.6g (%4.2f) %12.6g (%.3fs)%n",
+ name, b,
rms1, scale * rms1, max1, scale * max1, e1.getMean() / EPS,
rms2, scale * rms2, max2, scale * max2, e2.getMean() / EPS,
1e-9 * time);
