Repository navigation
fix overflow and precision loss in AccurateMath.IEEEremainder - #349
Conversation
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.
| // 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); |
There was a problem hiding this comment.
Can this use 2*p instead of p+p
There was a problem hiding this comment.
Yes, it's the same value. Doubling a double is exact, and 2 * p overflows to infinity for exactly the same p as p + p, so the overflow path where % leaves x unchanged behaves as before. Switched in 319201c and AccurateMathTest still passes. The two x + x comparisons below are the same case if you'd prefer those as 2 * x as well.
There was a problem hiding this comment.
Please change all cases of (x + x) to 2 * x. Thanks.
There was a problem hiding this comment.
Done in 4ddfabc. Both comparisons and the comment above them now use 2 * x, and AccurateMathTest still passes.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #349 +/- ##
============================================
+ Coverage 86.54% 87.22% +0.68%
+ Complexity 9787 91 -9696
============================================
Files 532 499 -33
Lines 35516 33538 -1978
Branches 6194 5854 -340
============================================
- Hits 30738 29255 -1483
+ Misses 3518 3161 -357
+ Partials 1260 1122 -138 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
With the default
JdkMathbackend (AccurateMath), on master:The general case evaluates
dividend - divisor * rint(dividend / divisor)in doubles. Past 2^53 the rounded quotient is not the nearest integer,divisor * nis itself rounded so the error scales with the dividend rather than the remainder, and a large quotient overflows the product. Even with |x/y| <= 1000, 974017 of 1000000 random pairs differ fromStrictMath.DerivativeStructure,SparseGradientandDSCompilerinherit this throughJdkMath.The fix reduces |dividend| modulo 2|divisor| with
%, which is exact for doubles, then subtracts |divisor| at most twice to round the quotient to even, as fdlibm's__ieee754_remainderdoes. Every step is exact and the method stays pure Java. It matchesStrictMath.IEEEremainderbit for bit on 100M pairs covering random bit patterns, subnormal divisors, ties and the overflow edges, on JDK 21 and 27. The javadoc bullet claiming the result always has the dividend's sign now applies to a zero result only, since IEEEremainder(3, 2) is -1.testIEEEremaindercompares against an exact BigDecimal reference; on master:mvn clean verifyis green across the reactor on JDK 21.