Match CPython remainder for subnormal inputs

src/m.rs

@@ -192,11 +192,6 @@ pub fn nextafter(x: f64, y: f64) -> f64 {
     unsafe { m_sys::nextafter(x, y) }
 }
 
-#[inline(always)]
-pub fn remainder(x: f64, y: f64) -> f64 {
-    unsafe { m_sys::remainder(x, y) }
-}
-
 #[inline(always)]
 pub fn trunc(n: f64) -> f64 {
     unsafe { m_sys::trunc(n) }

src/m_rust.rs

@@ -190,11 +190,6 @@ pub fn nextafter(x: f64, y: f64) -> f64 {
     libm::nextafter(x, y)
 }
 
-#[inline(always)]
-pub fn remainder(x: f64, y: f64) -> f64 {
-    libm::remainder(x, y)
-}
-
 #[inline(always)]
 pub fn trunc(n: f64) -> f64 {
     libm::trunc(n)

src/m_sys.rs

@@ -55,7 +55,6 @@ unsafe extern "C" {
     pub fn ldexp(x: f64, n: i32) -> f64;
     pub fn modf(n: f64, iptr: *mut f64) -> f64;
     pub fn nextafter(x: f64, y: f64) -> f64;
-    pub fn remainder(x: f64, y: f64) -> f64;
     pub fn trunc(n: f64) -> f64;
 
     // Special functions

src/math/misc.rs

@@ -247,12 +247,46 @@ pub fn fmod(x: f64, y: f64) -> Result<f64> {
 ///
 /// CPython implements this from scratch using fmod (m_remainder in
 /// mathmodule.c) rather than calling the C library's remainder().
-/// We delegate to libm's remainder() which is correct on all platforms
-/// where it conforms to IEEE 754. If you find a platform where the
-/// results differ from CPython, please file a bug.
 #[inline]
 pub fn remainder(x: f64, y: f64) -> Result<f64> {
-    super::math_2(x, y, crate::m::remainder)
+    super::math_2(x, y, m_remainder)
+}
+
+fn m_remainder(x: f64, y: f64) -> f64 {
+    if x.is_finite() && y.is_finite() {
+        if y == 0.0 {
+            return f64::NAN;
+        }
+
+        let absx = x.abs();
+        let absy = y.abs();
+        let m = m::fmod(absx, absy);
+        // Comparing with absy / 2 can round or underflow for subnormals.
+        // CPython compares the remainder with its complement instead.
+        let c = absy - m;
+        let r = if m < c {
+            m
+        } else if m > c {
+            -c
+        } else {
+            // At a tie, reducing half the preceding multiple modulo absy
+            // distinguishes odd and even multiples without rounding.
+            m - 2.0 * m::fmod(0.5 * (absx - m), absy)
+        };
+        return m::copysign(1.0, x) * r;
+    }
+
+    if x.is_nan() {
+        return x;
+    }
+    if y.is_nan() {
+        return y;
+    }
+    if x.is_infinite() {
+        return f64::NAN;
+    }
+    debug_assert!(y.is_infinite());
+    x
 }
 
 /// Return the value of the least significant bit of x.
@@ -474,6 +508,32 @@ mod tests {
         }
     }
 
+    #[test]
+    fn regression_remainder_subnormal_cases() {
+        let tiny = f64::from_bits(1);
+        for (x, y, expected) in [
+            (12.0, 23.0, -11.0),
+            (11.0, 23.0, 11.0),
+            (35.0, 23.0, -11.0),
+            (3.0, 2.0, -1.0),
+            (5.0, 2.0, 1.0),
+            (4.0, 2.0, 0.0),
+        ] {
+            for sx in [1.0, -1.0] {
+                for sy in [1.0, -1.0] {
+                    let x = sx * x * tiny;
+                    let y = sy * y * tiny;
+                    assert_eq!(
+                        remainder(x, y).unwrap().to_bits(),
+                        (sx * expected * tiny).to_bits(),
+                        "remainder({x}, {y})",
+                    );
+                    test_remainder(x, y);
+                }
+            }
+        }
+    }
+
     #[test]
     fn regression_remainder_halfway_even_cases() {
         // These cases exercise the half-way branch in CPython's custom
@@ -574,6 +634,18 @@ mod tests {
             test_remainder(x, y);
         }
 
+        #[test]
+        fn proptest_remainder_subnormal(
+            x in 1_u64..(1_u64 << 52),
+            y in 1_u64..(1_u64 << 52),
+            negative_x: bool,
+            negative_y: bool,
+        ) {
+            let x = f64::from_bits(x | ((negative_x as u64) << 63));
+            let y = f64::from_bits(y | ((negative_y as u64) << 63));
+            test_remainder(x, y);
+        }
+
         #[test]
         fn proptest_copysign(x: f64, y: f64) {
             test_copysign(x, y);