fix overflow and precision loss in AccurateMath.IEEEremainder

commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java

@@ -3691,7 +3691,7 @@ public static double hypot(final double x, final double y) {
      * {@code n} is the integer that is even.
      * <ul>
      * <li>If either operand is NaN, the result is NaN.</li>
-     * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
+     * <li>If the result is zero, its sign equals the sign of the dividend.</li>
      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
@@ -3710,9 +3710,20 @@ public static double IEEEremainder(final double dividend, final double divisor)
             }
         } else {
             // we are in the general case
-            final double n         = AccurateMath.rint(dividend / divisor);
-            final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
-            return (remainder == 0) ? AccurateMath.copySign(remainder, dividend) : remainder;
+            // Reduce |dividend| modulo 2 |divisor| with the exact truncating remainder
+            // (x is unchanged if 2 |divisor| overflows), then subtract |divisor| at most
+            // twice to round the quotient to even, as in fdlibm __ieee754_remainder.
+            // All steps are exact: x + x only overflows when x + x > p holds anyway.
+            final double p = Math.abs(divisor);
+            double x = Math.abs(dividend) % (2 * p);
+            if (x + x > p) {
+                x -= p;
+                if (x + x >= p) {
+                    x -= p;
+                }
+            }
+            // Restore the sign of the dividend (including -0.0)
+            return Double.doubleToRawLongBits(dividend) < 0 ? -x : x;
         }
     }
 

commons-math-legacy-core/src/test/java/org/apache/commons/math4/legacy/core/jdkmath/AccurateMathTest.java

@@ -1358,6 +1358,59 @@ public void testLogWithBase() {
         assertEquals(Double.NEGATIVE_INFINITY, AccurateMath.log(10, 0), 0);
     }
 
+    @Test
+    public void testIEEEremainder() {
+        // Special cases taken from the javadoc
+        Assertions.assertEquals(Double.NaN, AccurateMath.IEEEremainder(Double.NaN, 1.0));
+        Assertions.assertEquals(Double.NaN, AccurateMath.IEEEremainder(1.0, Double.NaN));
+        Assertions.assertEquals(Double.NaN, AccurateMath.IEEEremainder(Double.POSITIVE_INFINITY, 1.0));
+        Assertions.assertEquals(Double.NaN, AccurateMath.IEEEremainder(1.0, -0.0));
+        Assertions.assertEquals(-1.0, AccurateMath.IEEEremainder(-1.0, Double.POSITIVE_INFINITY));
+        Assertions.assertEquals(-0.0, AccurateMath.IEEEremainder(-0.0, 1.0));
+        Assertions.assertEquals(-0.0, AccurateMath.IEEEremainder(-4.0, 2.0));
+        // Quotient rounded to the nearest even integer
+        Assertions.assertEquals(-1.0, AccurateMath.IEEEremainder(3.0, 2.0));
+        Assertions.assertEquals(1.0, AccurateMath.IEEEremainder(5.0, 2.0));
+
+        // x - y * rint(x / y) overflows, or loses the remainder when y * n is inexact
+        Assertions.assertEquals(0.5, AccurateMath.IEEEremainder(Double.MAX_VALUE, 1.5));
+        Assertions.assertEquals(1.0, AccurateMath.IEEEremainder(-Double.MAX_VALUE, 3.0));
+        Assertions.assertEquals(-0.5, AccurateMath.IEEEremainder(0x1p52, 1.5));
+        Assertions.assertEquals(1.6653345369377348E-16, AccurateMath.IEEEremainder(1.5, 0.49999999999999994));
+
+        for (int i = 0; i < NUMBER_OF_TRIALS; i++) {
+            final double y = randomFiniteDouble();
+            // Alternate between unrestricted quotients and quotients below 1000
+            final double x = (i & 1) == 0 ? randomFiniteDouble() : y * (generator.nextDouble() * 2000 - 1000);
+            if (y == 0 || !Double.isFinite(x)) {
+                continue;
+            }
+            Assertions.assertEquals(exactIEEEremainder(x, y), AccurateMath.IEEEremainder(x, y),
+                () -> "IEEEremainder(" + x + ", " + y + ")");
+        }
+    }
+
+    private double randomFiniteDouble() {
+        double x;
+        do {
+            x = Double.longBitsToDouble(generator.nextLong());
+        } while (!Double.isFinite(x));
+        return x;
+    }
+
+    /**
+     * Reference IEEE remainder {@code x - y * n} with {@code n} the integer nearest
+     * to {@code x / y} (ties to even), computed exactly.
+     */
+    private static double exactIEEEremainder(double x, double y) {
+        final BigDecimal bx = new BigDecimal(x);
+        final BigDecimal by = new BigDecimal(y);
+        final BigDecimal n = bx.divide(by, 0, RoundingMode.HALF_EVEN);
+        // The remainder is exactly representable as a double
+        final double r = bx.subtract(by.multiply(n)).doubleValue();
+        return r == 0 ? Math.copySign(0.0, x) : r;
+    }
+
     @Test
     public void testIndicatorDouble() {
         double delta = 0.0;