Update remainder module of math

- Implement remainder function with test case
  - math.remainder was added to CPython in 3.7 and RustPython CI runs on 3.6

tests/snippets/stdlib_math.py

@@ -265,4 +265,95 @@ assert math.fmod(-3.0, INF) == -3.0
 assert math.fmod(3.0, NINF) == 3.0
 assert math.fmod(-3.0, NINF) == -3.0
 assert math.fmod(0.0, 3.0) == 0.0
-assert math.fmod(0.0, NINF) == 0.0
+assert math.fmod(0.0, NINF) == 0.0
+
+"""
+TODO: math.remainder was added to CPython in 3.7 and RustPython CI runs on 3.6.
+So put the tests of math.remainder in a comment for now.
+https://github.com/RustPython/RustPython/pull/1589#issuecomment-551424940
+"""
+
+# testcases = [
+#     # Remainders modulo 1, showing the ties-to-even behaviour.
+#     '-4.0 1 -0.0',
+#     '-3.8 1  0.8',
+#     '-3.0 1 -0.0',
+#     '-2.8 1 -0.8',
+#     '-2.0 1 -0.0',
+#     '-1.8 1  0.8',
+#     '-1.0 1 -0.0',
+#     '-0.8 1 -0.8',
+#     '-0.0 1 -0.0',
+#     ' 0.0 1  0.0',
+#     ' 0.8 1  0.8',
+#     ' 1.0 1  0.0',
+#     ' 1.8 1 -0.8',
+#     ' 2.0 1  0.0',
+#     ' 2.8 1  0.8',
+#     ' 3.0 1  0.0',
+#     ' 3.8 1 -0.8',
+#     ' 4.0 1  0.0',
+
+#     # Reductions modulo 2*pi
+#     '0x0.0p+0 0x1.921fb54442d18p+2 0x0.0p+0',
+#     '0x1.921fb54442d18p+0 0x1.921fb54442d18p+2  0x1.921fb54442d18p+0',
+#     '0x1.921fb54442d17p+1 0x1.921fb54442d18p+2  0x1.921fb54442d17p+1',
+#     '0x1.921fb54442d18p+1 0x1.921fb54442d18p+2  0x1.921fb54442d18p+1',
+#     '0x1.921fb54442d19p+1 0x1.921fb54442d18p+2 -0x1.921fb54442d17p+1',
+#     '0x1.921fb54442d17p+2 0x1.921fb54442d18p+2 -0x0.0000000000001p+2',
+#     '0x1.921fb54442d18p+2 0x1.921fb54442d18p+2  0x0p0',
+#     '0x1.921fb54442d19p+2 0x1.921fb54442d18p+2  0x0.0000000000001p+2',
+#     '0x1.2d97c7f3321d1p+3 0x1.921fb54442d18p+2  0x1.921fb54442d14p+1',
+#     '0x1.2d97c7f3321d2p+3 0x1.921fb54442d18p+2 -0x1.921fb54442d18p+1',
+#     '0x1.2d97c7f3321d3p+3 0x1.921fb54442d18p+2 -0x1.921fb54442d14p+1',
+#     '0x1.921fb54442d17p+3 0x1.921fb54442d18p+2 -0x0.0000000000001p+3',
+#     '0x1.921fb54442d18p+3 0x1.921fb54442d18p+2  0x0p0',
+#     '0x1.921fb54442d19p+3 0x1.921fb54442d18p+2  0x0.0000000000001p+3',
+#     '0x1.f6a7a2955385dp+3 0x1.921fb54442d18p+2  0x1.921fb54442d14p+1',
+#     '0x1.f6a7a2955385ep+3 0x1.921fb54442d18p+2  0x1.921fb54442d18p+1',
+#     '0x1.f6a7a2955385fp+3 0x1.921fb54442d18p+2 -0x1.921fb54442d14p+1',
+#     '0x1.1475cc9eedf00p+5 0x1.921fb54442d18p+2  0x1.921fb54442d10p+1',
+#     '0x1.1475cc9eedf01p+5 0x1.921fb54442d18p+2 -0x1.921fb54442d10p+1',
+
+#     # Symmetry with respect to signs.
+#     ' 1  0.c  0.4',
+#     '-1  0.c -0.4',
+#     ' 1 -0.c  0.4',
+#     '-1 -0.c -0.4',
+#     ' 1.4  0.c -0.4',
+#     '-1.4  0.c  0.4',
+#     ' 1.4 -0.c -0.4',
+#     '-1.4 -0.c  0.4',
+
+#     # Huge modulus, to check that the underlying algorithm doesn't
+#     # rely on 2.0 * modulus being representable.
+#     '0x1.dp+1023 0x1.4p+1023  0x0.9p+1023',
+#     '0x1.ep+1023 0x1.4p+1023 -0x0.ap+1023',
+#     '0x1.fp+1023 0x1.4p+1023 -0x0.9p+1023',
+# ]
+
+# for case in testcases:
+#     x_hex, y_hex, expected_hex = case.split()
+#     # print(x_hex, y_hex, expected_hex)
+#     x = float.fromhex(x_hex)
+#     y = float.fromhex(y_hex)
+#     expected = float.fromhex(expected_hex)
+#     actual = math.remainder(x, y)
+#     # Cheap way of checking that the floats are
+#     # as identical as we need them to be.
+#     assert actual.hex() == expected.hex()
+#     # self.assertEqual(actual.hex(), expected.hex())
+
+
+# # Test tiny subnormal modulus: there's potential for
+# # getting the implementation wrong here (for example,
+# # by assuming that modulus/2 is exactly representable).
+# tiny = float.fromhex('1p-1074')  # min +ve subnormal
+# for n in range(-25, 25):
+#     if n == 0:
+#         continue
+#     y = n * tiny
+#     for m in range(100):
+#         x = m * tiny
+#         actual = math.remainder(x, y)
+#         actual = math.remainder(-x, y)

vm/src/stdlib/math.rs

@@ -275,13 +275,20 @@ fn math_modf(x: IntoPyFloat, _vm: &VirtualMachine) -> (f64, f64) {
     (x.fract(), x.trunc())
 }
 
+fn fmod(x: f64, y: f64) -> f64 {
+    if y.is_infinite() && x.is_finite() {
+        return x;
+    }
+
+    x % y
+}
+
 fn math_fmod(x: IntoPyFloat, y: IntoPyFloat, vm: &VirtualMachine) -> PyResult<f64> {
     let x = x.to_f64();
     let y = y.to_f64();
-    if y.is_infinite() && x.is_finite() {
-        return Ok(x);
-    }
-    let r = x % y;
+
+    let r = fmod(x, y);
+
     if r.is_nan() && !x.is_nan() && !y.is_nan() {
         return Err(vm.new_value_error("math domain error".to_string()));
     }
@@ -289,6 +296,46 @@ fn math_fmod(x: IntoPyFloat, y: IntoPyFloat, vm: &VirtualMachine) -> PyResult<f6
     Ok(r)
 }
 
+fn math_remainder(x: IntoPyFloat, y: IntoPyFloat, vm: &VirtualMachine) -> PyResult<f64> {
+    let x = x.to_f64();
+    let y = y.to_f64();
+    if x.is_finite() && y.is_finite() {
+        if y == 0.0 {
+            return Ok(std::f64::NAN);
+        }
+
+        let absx = x.abs();
+        let absy = y.abs();
+        let modulus = absx % absy;
+
+        let c = absy - modulus;
+        let r;
+        if modulus < c {
+            r = modulus;
+        } else if modulus > c {
+            r = -c;
+        } else {
+            r = modulus - 2.0 * fmod(0.5 * (absx - modulus), absy);
+        }
+
+        return Ok(1.0_f64.copysign(x) * r);
+    }
+
+    if x.is_nan() {
+        return Ok(x);
+    }
+    if y.is_nan() {
+        return Ok(y);
+    }
+    if x.is_infinite() {
+        return Ok(std::f64::NAN);
+    }
+    if y.is_infinite() {
+        return Err(vm.new_value_error("math domain error".to_string()));
+    }
+    Ok(x)
+}
+
 pub fn make_module(vm: &VirtualMachine) -> PyObjectRef {
     let ctx = &vm.ctx;
 
@@ -342,6 +389,7 @@ pub fn make_module(vm: &VirtualMachine) -> PyObjectRef {
         "ldexp" => ctx.new_rustfunc(math_ldexp),
         "modf" => ctx.new_rustfunc(math_modf),
         "fmod" => ctx.new_rustfunc(math_fmod),
+        "remainder" => ctx.new_rustfunc(math_remainder),
 
         // Rounding functions:
         "trunc" => ctx.new_rustfunc(math_trunc),