add `cpy11__float_div_mod`

docs/modules/math.md

@@ -123,11 +123,14 @@ Convert angle `x` from radians to degrees.
 
 Convert angle `x` from degrees to radians.
 
-
 ### `math.modf(x)`
 
 Return the fractional and integer parts of `x`. Both results carry the sign of `x` and are floats.
 
+### `math.copysign(x, y)`
+
+Return a float with the magnitude (absolute value) of `x` but the sign of `y`.
+
 ### `math.factorial(x)`
 
-Return `x` factorial as an integer.
+Return `x` factorial as an integer.

src/bindings/py_number.c

@@ -149,6 +149,44 @@ static py_i64 cpy11__fast_mod(py_i64 a, py_i64 b) {
     return b < 0 ? -res : res;
 }
 
+// https://github.com/python/cpython/blob/3.11/Objects/floatobject.c#L677
+static void cpy11__float_div_mod(double vx, double wx, double *floordiv, double *mod)
+{
+    double div;
+    *mod = fmod(vx, wx);
+    /* fmod is typically exact, so vx-mod is *mathematically* an
+       exact multiple of wx.  But this is fp arithmetic, and fp
+       vx - mod is an approximation; the result is that div may
+       not be an exact integral value after the division, although
+       it will always be very close to one.
+    */
+    div = (vx - *mod) / wx;
+    if (*mod) {
+        /* ensure the remainder has the same sign as the denominator */
+        if ((wx < 0) != (*mod < 0)) {
+            *mod += wx;
+            div -= 1.0;
+        }
+    }
+    else {
+        /* the remainder is zero, and in the presence of signed zeroes
+           fmod returns different results across platforms; ensure
+           it has the same sign as the denominator. */
+        *mod = copysign(0.0, wx);
+    }
+    /* snap quotient to nearest integral value */
+    if (div) {
+        *floordiv = floor(div);
+        if (div - *floordiv > 0.5) {
+            *floordiv += 1.0;
+        }
+    }
+    else {
+        /* div is zero - get the same sign as the true quotient */
+        *floordiv = copysign(0.0, vx / wx); /* zero w/ sign of vx/wx */
+    }
+}
+
 static bool int__floordiv__(int argc, py_Ref argv) {
     PY_CHECK_ARGC(2);
     py_i64 lhs = py_toint(&argv[0]);
@@ -181,8 +219,24 @@ static bool float__floordiv__(int argc, py_Ref argv) {
     py_f64 rhs;
     if(try_castfloat(&argv[1], &rhs)) {
         if(rhs == 0.0) return ZeroDivisionError("float modulo by zero");
-        py_f64 r = fmod(lhs, rhs);
-        py_newfloat(py_retval(), trunc((lhs - r) / rhs));
+        double q, r;
+        cpy11__float_div_mod(lhs, rhs, &q, &r);
+        py_newfloat(py_retval(), q);
+        return true;
+    }
+    py_newnotimplemented(py_retval());
+    return true;
+}
+
+static bool float__rfloordiv__(int argc, py_Ref argv) {
+    PY_CHECK_ARGC(2);
+    py_f64 rhs = py_tofloat(&argv[0]);
+    py_f64 lhs;
+    if(try_castfloat(&argv[1], &lhs)) {
+        if(rhs == 0.0) return ZeroDivisionError("float modulo by zero");
+        double q, r;
+        cpy11__float_div_mod(lhs, rhs, &q, &r);
+        py_newfloat(py_retval(), q);
         return true;
     }
     py_newnotimplemented(py_retval());
@@ -195,7 +249,9 @@ static bool float__mod__(int argc, py_Ref argv) {
     py_f64 rhs;
     if(try_castfloat(&argv[1], &rhs)) {
         if(rhs == 0.0) return ZeroDivisionError("float modulo by zero");
-        py_newfloat(py_retval(), fmod(lhs, rhs));
+        double q, r;
+        cpy11__float_div_mod(lhs, rhs, &q, &r);
+        py_newfloat(py_retval(), r);
         return true;
     }
     py_newnotimplemented(py_retval());
@@ -208,7 +264,9 @@ static bool float__rmod__(int argc, py_Ref argv) {
     py_f64 lhs;
     if(try_castfloat(&argv[1], &lhs)) {
         if(rhs == 0.0) return ZeroDivisionError("float modulo by zero");
-        py_newfloat(py_retval(), fmod(lhs, rhs));
+        double q, r;
+        cpy11__float_div_mod(lhs, rhs, &q, &r);
+        py_newfloat(py_retval(), r);
         return true;
     }
     py_newnotimplemented(py_retval());
@@ -221,9 +279,10 @@ static bool float__divmod__(int argc, py_Ref argv) {
     py_f64 rhs;
     if(try_castfloat(&argv[1], &rhs)) {
         if(rhs == 0.0) return ZeroDivisionError("float modulo by zero");
-        py_f64 r = fmod(lhs, rhs);
+        double q, r;
+        cpy11__float_div_mod(lhs, rhs, &q, &r);
         py_Ref p = py_newtuple(py_retval(), 2);
-        py_newfloat(&p[0], trunc((lhs - r) / rhs));
+        py_newfloat(&p[0], q);
         py_newfloat(&p[1], r);
         return true;
     }
@@ -565,6 +624,7 @@ void pk_number__register() {
 
     // fmod
     py_bindmagic(tp_float, __floordiv__, float__floordiv__);
+    py_bindmagic(tp_float, __rfloordiv__, float__rfloordiv__);
     py_bindmagic(tp_float, __mod__, float__mod__);
     py_bindmagic(tp_float, __rmod__, float__rmod__);
     py_bindmagic(tp_float, __divmod__, float__divmod__);

src/modules/math.c

@@ -133,6 +133,7 @@ static bool math_radians(int argc, py_Ref argv) {
 }
 
 TWO_ARG_FUNC(fmod, fmod)
+TWO_ARG_FUNC(copysign, copysign)
 
 static bool math_modf(int argc, py_Ref argv) {
     PY_CHECK_ARGC(1);
@@ -200,6 +201,7 @@ void pk__add_module_math() {
 
     py_bindfunc(mod, "fmod", math_fmod);
     py_bindfunc(mod, "modf", math_modf);
+    py_bindfunc(mod, "copysign", math_copysign);
     py_bindfunc(mod, "factorial", math_factorial);
 }
 

tests/02_float.py

@@ -120,3 +120,10 @@ assert eq(10.5 // 4.5, 2.0)
 _0, _1 = divmod(10.5, 4)
 assert eq(_0, 2.0)
 assert eq(_1, 2.5)
+
+assert eq(3.4 % -2, -0.6)
+assert eq(-2 % 3.4, 1.4)
+assert eq(-3.4 % -2, -1.4)
+assert eq(-6 // 3.4, -2.0)
+assert eq(-6 % 3.4, 0.8)
+