ENH: Make numpy ufuncs compatible with Python divmod.

The following numpy ufuncs are reimplemented using the npy_divmod
function.

- remainder ('%')
- floor_divide ('//')

Example of differences

Currently

    In [1]: a = np.array(78 * 6e-8)

    In [2]: b = np.array(6e-8)

    In [3]: a // b
    Out[3]: 77.0

    In [4]: a % b
    Out[4]: 5.9999999999999651e-08

Previously

    In [1]: a = np.array(78 * 6e-8)

    In [2]: b = np.array(6e-8)

    In [3]: a // b
    Out[3]: 78.0

    In [4]: a % b
    Out[4]: 0.0

The first agrees with the Python values for the same operation and is a
bit more accurate for the input values as represented internally.

Closes #7224.

1 files changed, 10 insertions, 15 deletions

diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 7b8dcdbaf..0d9806f5d 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -1696,7 +1696,8 @@ NPY_NO_EXPORT void
     BINARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
         const @type@ in2 = *(@type@ *)ip2;
-        *((@type@ *)op1) = npy_floor@c@(in1/in2);
+        @type@ mod;
+        *((@type@ *)op1) = npy_divmod@c@(in1, in2, &mod);
     }
 }
 
@@ -1706,8 +1707,7 @@ NPY_NO_EXPORT void
     BINARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
         const @type@ in2 = *(@type@ *)ip2;
-        const @type@ div = in1/in2;
-        *((@type@ *)op1) = in2*(div - npy_floor@c@(div));
+        npy_divmod@c@(in1, in2, (@type@ *)op1);
     }
 }
 
@@ -2011,9 +2011,10 @@ NPY_NO_EXPORT void
 HALF_floor_divide(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     BINARY_LOOP {
-        const float in1 = npy_half_to_float(*(npy_half *)ip1);
-        const float in2 = npy_half_to_float(*(npy_half *)ip2);
-        *((npy_half *)op1) = npy_float_to_half(npy_floorf(in1/in2));
+        const npy_half in1 = *(npy_half *)ip1;
+        const npy_half in2 = *(npy_half *)ip2;
+        npy_half mod;
+        *((npy_half *)op1) = npy_half_divmod(in1, in2, &mod);
     }
 }
 
@@ -2021,15 +2022,9 @@ NPY_NO_EXPORT void
 HALF_remainder(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     BINARY_LOOP {
-        const float in1 = npy_half_to_float(*(npy_half *)ip1);
-        const float in2 = npy_half_to_float(*(npy_half *)ip2);
-        const float res = npy_fmodf(in1,in2);
-        if (res && ((in2 < 0) != (res < 0))) {
-            *((npy_half *)op1) = npy_float_to_half(res + in2);
-        }
-        else {
-            *((npy_half *)op1) = npy_float_to_half(res);
-        }
+        const npy_half in1 = *(npy_half *)ip1;
+        const npy_half in2 = *(npy_half *)ip2;
+        npy_half_divmod(in1, in2, (npy_half *)op1);
     }
 }