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 2023/08/30 11:29:59 UTC
[commons-numbers] branch master updated: NUMBERS-204: Correct sub-normal round-off computation in sum of products
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 c8457580 NUMBERS-204: Correct sub-normal round-off computation in sum of products
c8457580 is described below
commit c8457580f049c7e8ad7507651ccd9bf8ff05aca5
Author: aherbert <ah...@apache.org>
AuthorDate: Wed Aug 30 12:25:03 2023 +0100
NUMBERS-204: Correct sub-normal round-off computation in sum of products
---
.../commons/numbers/core/ExtendedPrecision.java | 25 +++++-
.../commons/numbers/core/DoubleTestUtils.java | 13 ++++
.../numbers/core/ExtendedPrecisionTest.java | 90 ++++++++++++++++++++--
.../org/apache/commons/numbers/core/SumTest.java | 15 ++++
.../org/apache/commons/numbers/core/TestUtils.java | 2 +-
src/changes/changes.xml | 3 +
6 files changed, 137 insertions(+), 11 deletions(-)
diff --git a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/ExtendedPrecision.java b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/ExtendedPrecision.java
index df256a29..5417d70f 100644
--- a/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/ExtendedPrecision.java
+++ b/commons-numbers-core/src/main/java/org/apache/commons/numbers/core/ExtendedPrecision.java
@@ -30,15 +30,24 @@ final class ExtendedPrecision {
* 996 is the value obtained from {@code Math.getExponent(Double.MAX_VALUE / MULTIPLIER)}.
*/
private static final double SAFE_UPPER = 0x1.0p996;
+ /**
+ * The lower limit for a product {@code x * y} below which the round-off component may be
+ * sub-normal. This is set as 2^-1022 * 2^54.
+ */
+ private static final double SAFE_LOWER = 0x1.0p-968;
/** The scale to use when down-scaling during a split into a high part.
* This must be smaller than the inverse of the multiplier and a power of 2 for exact scaling. */
private static final double DOWN_SCALE = 0x1.0p-30;
-
- /** The scale to use when re-scaling during a split into a high part.
+ /** The scale to use when up-scaling during a split into a high part.
* This is the inverse of {@link #DOWN_SCALE}. */
private static final double UP_SCALE = 0x1.0p30;
+ /** The upscale factor squared. */
+ private static final double UP_SCALE2 = 0x1.0p60;
+ /** The downscale factor squared. */
+ private static final double DOWN_SCALE2 = 0x1.0p-60;
+
/** Private constructor. */
private ExtendedPrecision() {
// intentionally empty.
@@ -103,7 +112,8 @@ final class ExtendedPrecision {
// 3-fold higher than any component.
final double a = Math.abs(x);
final double b = Math.abs(y);
- if (a + b + Math.abs(xy) >= SAFE_UPPER) {
+ final double ab = Math.abs(xy);
+ if (a + b + ab >= SAFE_UPPER) {
// Only required to scale the largest number as x*y does not overflow.
if (a > b) {
return DD.twoProductLow(x * DOWN_SCALE, y, xy * DOWN_SCALE) * UP_SCALE;
@@ -111,6 +121,15 @@ final class ExtendedPrecision {
return DD.twoProductLow(x, y * DOWN_SCALE, xy * DOWN_SCALE) * UP_SCALE;
}
+ // The result is computed using a product of the low parts.
+ // To avoid underflow in the low parts we note that these are approximately a factor
+ // of 2^27 smaller than the original inputs so their product will be ~2^54 smaller
+ // than the product xy. Ensure the product is at least 2^54 above a sub-normal.
+ if (ab <= SAFE_LOWER) {
+ // Scaling up here is safe: the largest magnitude cannot be above SAFE_LOWER / MIN_VALUE.
+ return DD.twoProductLow(x * UP_SCALE, y * UP_SCALE, xy * UP_SCALE2) * DOWN_SCALE2;
+ }
+
// No scaling required
return DD.twoProductLow(x, y, xy);
}
diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DoubleTestUtils.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DoubleTestUtils.java
index 5546ab70..4cf9b99c 100644
--- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DoubleTestUtils.java
+++ b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/DoubleTestUtils.java
@@ -66,4 +66,17 @@ final class DoubleTestUtils {
final long exp = rng.nextInt(expRange) + minExp + 1023;
return Double.longBitsToDouble(bits | (exp << 52));
}
+
+ /** Construct a random double with an exponent of 0.
+ * The magnitude will be in the range {@code [1, 2)} and the sign is random.
+ * @param rng random number generator
+ * @return random double in +/-[1, 2)
+ */
+ public static double randomDouble(final UniformRandomProvider rng) {
+ // Create random doubles using random bits in the sign bit and the mantissa.
+ final long mask = ((1L << 52) - 1) | 1L << 63;
+ final long bits = rng.nextLong() & mask;
+ // The exponent must be unsigned so + 1023 to the signed exponent
+ return Double.longBitsToDouble(bits | (1023L << 52));
+ }
}
diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/ExtendedPrecisionTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/ExtendedPrecisionTest.java
index 0cfe0753..1c19e3bf 100644
--- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/ExtendedPrecisionTest.java
+++ b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/ExtendedPrecisionTest.java
@@ -16,8 +16,15 @@
*/
package org.apache.commons.numbers.core;
+import java.util.function.Supplier;
+import java.util.stream.Stream;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.Arguments;
+import org.junit.jupiter.params.provider.MethodSource;
/**
* Test cases for the {@link ExtendedPrecision} class.
@@ -51,10 +58,14 @@ class ExtendedPrecisionTest {
/**
* Test {@link ExtendedPrecision#productLow(double, double, double)} computes the same
- * result as JDK 9 Math.fma(x, y, -x * y) for edge cases.
+ * result as JDK 9 Math.fma(x, y, -x * y) for edge cases, e.g.
+ * <pre>
+ * jshell> static double f(double a, double b) { return Math.fma(x, y, -x * y); }
+ * jshell> f(5.6266120027810604E-148, 2.9150607442566245E-154);
+ * </pre>
*/
@Test
- void testProductLow() {
+ void testProductLowEdgeCases() {
assertProductLow(0.0, 1.0, Math.nextDown(Double.MIN_NORMAL));
assertProductLow(0.0, -1.0, Math.nextDown(Double.MIN_NORMAL));
assertProductLow(Double.NaN, 1.0, Double.POSITIVE_INFINITY);
@@ -62,11 +73,10 @@ class ExtendedPrecisionTest {
assertProductLow(Double.NaN, 1.0, Double.NaN);
assertProductLow(0.0, 1.0, Double.MAX_VALUE);
assertProductLow(Double.NaN, 2.0, Double.MAX_VALUE);
- }
-
- private static void assertProductLow(double expected, double x, double y) {
- // Requires a delta of 0.0 to assert -0.0 == 0.0
- Assertions.assertEquals(expected, ExtendedPrecision.productLow(x, y, x * y), 0.0);
+ // Product is normal, round-off is sub-normal
+ // Ignoring sub-normals during computation results in a rounding error
+ assertProductLow(-0.0, -2.73551683292218E-154, -1.0861547555023299E-154);
+ assertProductLow(9.023244E-318, 5.6266120027810604E-148, 2.9150607442566245E-154);
}
/**
@@ -92,4 +102,70 @@ class ExtendedPrecisionTest {
Assertions.assertTrue(Math.abs(hi2) > Math.abs(lo2));
}
+
+ @ParameterizedTest
+ @MethodSource
+ void testProductLow(double x, double y) {
+ // Assumes the arguments are in [1, 2) so scaling factors hit target cases
+ Assertions.assertTrue(Math.abs(x) >= 1 && Math.abs(x) < 2, () -> "Invalid x: " + x);
+ Assertions.assertTrue(Math.abs(y) >= 1 && Math.abs(y) < 2, () -> "Invalid y: " + y);
+
+ final double low = DD.ofProduct(x, y).lo();
+ assertProductLow(low, x, y);
+
+ // Product approaching and into sub-normal
+ final int[] scaleDown = {-490, -510, -511, -512, -513};
+ for (final int e1 : scaleDown) {
+ final double a = Math.scalb(x, e1);
+ for (final int e2 : scaleDown) {
+ final double b = Math.scalb(y, e2);
+ final double expected = Math.scalb(low, e1 + e2);
+ assertProductLow(expected, a, b, () -> "Product towards sub-normal");
+ }
+ }
+
+ // Product approaching overflow
+ final int[] scaleUp = {509, 510, 511, 512};
+ for (final int e1 : scaleUp) {
+ final double a = Math.scalb(x, e1);
+ for (final int e2 : scaleUp) {
+ final double b = Math.scalb(y, e2);
+ if (Double.isFinite(a * b)) {
+ final double expected = Math.scalb(low, e1 + e2);
+ assertProductLow(expected, a, b, () -> "Product towards overflow");
+ }
+ }
+ }
+
+ // Split of an argument approaching overflow
+ final int[] scaleUp2 = {990, 1000, 1010};
+ for (final int e : scaleUp2) {
+ final double a = Math.scalb(x, e);
+ final double expected = Math.scalb(low, e);
+ assertProductLow(expected, a, y, () -> "Argument x split towards overflow");
+ assertProductLow(expected, y, a, () -> "Argument y split towards overflow");
+ }
+ }
+
+ static Stream<Arguments> testProductLow() {
+ final Stream.Builder<Arguments> builder = Stream.builder();
+ final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
+ final int n = 50;
+ // Generate arguments in [1, 2)
+ for (int i = 0; i < n; i++) {
+ builder.add(Arguments.of(DoubleTestUtils.randomDouble(rng),
+ DoubleTestUtils.randomDouble(rng)));
+ }
+ return builder.build();
+ }
+
+ private static void assertProductLow(double expected, double x, double y) {
+ assertProductLow(expected, x, y, null);
+ }
+
+ private static void assertProductLow(double expected, double x, double y, Supplier<String> msg) {
+ // Allowed ULP=0. This allows -0.0 and 0.0 to be equal.
+ TestUtils.assertEquals(expected, ExtendedPrecision.productLow(x, y, x * y), 0,
+ () -> TestUtils.prefix(msg) + "low(" + x + " * " + y + ")");
+ }
}
diff --git a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java
index 844d7096..e29fc858 100644
--- a/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java
+++ b/commons-numbers-core/src/test/java/org/apache/commons/numbers/core/SumTest.java
@@ -25,6 +25,7 @@ import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
import org.junit.jupiter.params.ParameterizedTest;
import org.junit.jupiter.params.provider.Arguments;
+import org.junit.jupiter.params.provider.CsvSource;
import org.junit.jupiter.params.provider.MethodSource;
class SumTest {
@@ -259,6 +260,20 @@ class SumTest {
Assertions.assertTrue(Math.abs(naive - sum) > 1.5);
}
+ @ParameterizedTest
+ @CsvSource({
+ // See: NUMBERS-204
+ // Round-off == 0.0
+ "-2.73551683292218E-154, -1.0861547555023299E-154",
+ // Round-off = 4.9E-324
+ "1.4134286753429383E-154, -4.1395762189346144E-154",
+ })
+ void testSumOfProduct_tiny(double x, double y) {
+ final Sum s = Sum.create();
+ s.addProduct(x, y);
+ Assertions.assertEquals(x * y, s.getAsDouble());
+ }
+
@Test
void testSumOfProducts_nonFinite() {
// arrange
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 2d629dcc..270771f4 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
@@ -357,7 +357,7 @@ public final class TestUtils {
* @param msg Message supplier
* @return the prefix
*/
- private static String prefix(Supplier<String> msg) {
+ static String prefix(Supplier<String> msg) {
return msg == null ? "" : msg.get() + ": ";
}
}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 59de54ad..6823ef34 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -56,6 +56,9 @@ If the output is not quite correct, check for invisible trailing spaces!
<release version="1.2" date="TBD" description="
New features, updates and bug fixes.
">
+ <action dev="aherbert" type="fix" issue="NUMBERS-204">
+ "Sum": Correct sub-normal round-off computation in sum of products.
+ </action>
<action dev="aherbert" type="add" issue="NUMBERS-193">
"DD": Add an extended precision floating-point number. A double-double (DD) number
is an unevaluated sum of two IEEE double precision numbers capable of representing