From cac545754e2d492fa3547e6a904b832f96f3c957 Mon Sep 17 00:00:00 2001
From: Mahammad Insaf <Insaf@digiscrypt.com>
Date: Tue, 3 Nov 2026 21:07:54 +0530
Subject: [PATCH 1/2] fix overflow and precision loss in
 AccurateMath.IEEEremainder

The general case computed dividend - divisor * rint(dividend / divisor) in
double arithmetic. Once the quotient passes 2^53 rint no longer sees the
nearest integer, divisor * n is itself rounded, and for large quotients the
product overflows, so IEEEremainder(Double.MAX_VALUE, 1.5) returned -Infinity.

Reduce |dividend| modulo 2 |divisor| with the exact % operator and subtract
|divisor| at most twice to round the quotient to even, as fdlibm does. The
result now matches StrictMath.IEEEremainder.
---
 .../math4/core/jdkmath/AccurateMath.java      | 19 +++++--
 .../legacy/core/jdkmath/AccurateMathTest.java | 53 +++++++++++++++++++
 2 files changed, 68 insertions(+), 4 deletions(-)

diff --git a/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java b/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java
index bcbff70ad..09eee2861 100644
--- a/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java
+++ b/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) % (p + 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;
         }
     }
 
diff --git a/commons-math-legacy-core/src/test/java/org/apache/commons/math4/legacy/core/jdkmath/AccurateMathTest.java b/commons-math-legacy-core/src/test/java/org/apache/commons/math4/legacy/core/jdkmath/AccurateMathTest.java
index d75f7693f..021b9e229 100644
--- a/commons-math-legacy-core/src/test/java/org/apache/commons/math4/legacy/core/jdkmath/AccurateMathTest.java
+++ b/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;

From 319201c036a42dc1c7e659747cf986cde678f8c4 Mon Sep 17 00:00:00 2001
From: Mahammad Insaf <Insaf@digiscrypt.com>
Date: Fri, 6 Nov 2026 04:52:11 +0530
Subject: [PATCH 2/2] AccurateMath.IEEEremainder: use 2 * p for the modulus

---
 .../org/apache/commons/math4/core/jdkmath/AccurateMath.java     | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java b/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java
index 09eee2861..13798f8a1 100644
--- a/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java
+++ b/commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java
@@ -3715,7 +3715,7 @@ public static double IEEEremainder(final double dividend, final double divisor)
             // 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) % (p + p);
+            double x = Math.abs(dividend) % (2 * p);
             if (x + x > p) {
                 x -= p;
                 if (x + x >= p) {