ENH: Add new npy_divmod function to npy_math.

The new function is taken from the Python version of float_divmod and
computes the result of floor_division and modulus together so that they
can be kept compatible. This should also result in the '//' and '%'
operators applied to np.float64 and Python float returning the same
values.

The intent is to implement the ufuncs floor_divide and remainder using
the npy_divmod so that their results will also match those of '//' and
'%' while providing consistency between the two.

Note that npy_divmod uses fmod, which is very slow. As a result, the
floor_division and remainder functions are about 4x slower than the
previous versions based on the floor function, but should be more
accurate.

4 files changed, 64 insertions, 0 deletions

diff --git a/numpy/core/include/numpy/halffloat.h b/numpy/core/include/numpy/halffloat.h
index 944f0ea34..ab0d221fb 100644
--- a/numpy/core/include/numpy/halffloat.h
+++ b/numpy/core/include/numpy/halffloat.h
@@ -37,6 +37,7 @@ int npy_half_signbit(npy_half h);
 npy_half npy_half_copysign(npy_half x, npy_half y);
 npy_half npy_half_spacing(npy_half h);
 npy_half npy_half_nextafter(npy_half x, npy_half y);
+npy_half npy_half_divmod(npy_half x, npy_half y, npy_half *modulus);
 
 /*
  * Half-precision constants
diff --git a/numpy/core/include/numpy/npy_math.h b/numpy/core/include/numpy/npy_math.h
index 3dae583f3..e76508de0 100644
--- a/numpy/core/include/numpy/npy_math.h
+++ b/numpy/core/include/numpy/npy_math.h
@@ -309,16 +309,20 @@ double npy_deg2rad(double x);
 double npy_rad2deg(double x);
 double npy_logaddexp(double x, double y);
 double npy_logaddexp2(double x, double y);
+double npy_divmod(double x, double y, double *modulus);
 
 float npy_deg2radf(float x);
 float npy_rad2degf(float x);
 float npy_logaddexpf(float x, float y);
 float npy_logaddexp2f(float x, float y);
+float npy_divmodf(float x, float y, float *modulus);
 
 npy_longdouble npy_deg2radl(npy_longdouble x);
 npy_longdouble npy_rad2degl(npy_longdouble x);
 npy_longdouble npy_logaddexpl(npy_longdouble x, npy_longdouble y);
 npy_longdouble npy_logaddexp2l(npy_longdouble x, npy_longdouble y);
+npy_longdouble npy_divmodl(npy_longdouble x, npy_longdouble y,
+                           npy_longdouble *modulus);
 
 #define npy_degrees npy_rad2deg
 #define npy_degreesf npy_rad2degf
diff --git a/numpy/core/src/npymath/halffloat.c b/numpy/core/src/npymath/halffloat.c
index 34ac64287..951768256 100644
--- a/numpy/core/src/npymath/halffloat.c
+++ b/numpy/core/src/npymath/halffloat.c
@@ -224,6 +224,17 @@ int npy_half_ge(npy_half h1, npy_half h2)
     return npy_half_le(h2, h1);
 }
 
+npy_half npy_half_divmod(npy_half h1, npy_half h2, npy_half *modulus)
+{
+    float fh1 = npy_half_to_float(h1);
+    float fh2 = npy_half_to_float(h2);
+    float div, mod;
+
+    div = npy_divmodf(fh1, fh2, &mod);
+    *modulus = npy_float_to_half(mod);
+    return npy_float_to_half(div);
+}
+
 
 
 /*
diff --git a/numpy/core/src/npymath/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src
index 4dcb01986..45b618a56 100644
--- a/numpy/core/src/npymath/npy_math.c.src
+++ b/numpy/core/src/npymath/npy_math.c.src
@@ -608,6 +608,54 @@ double npy_log2(double x)
     }
 }
 
+/*
+ * Python version of divmod.
+ *
+ * The implementation is mostly copied from cpython 3.5.
+ */
+@type@
+npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus)
+{
+    @type@ div, mod, floordiv;
+
+    mod = npy_fmod@c@(a, b);
+
+    if (!b) {
+        /* If b == 0, return result of fmod. For IEEE is nan */
+        *modulus = mod;
+        return mod;
+    }
+
+    /* a - mod should be very nearly an integer multiple of b */
+    div = (a - mod) / b;
+
+    /* adjust fmod result to conform to Python convention of remainder */
+    if (mod) {
+        if ((b < 0) != (mod < 0)) {
+            mod += b;
+            div -= 1.0@c@;
+        }
+    }
+    else {
+        /* if mod is zero ensure correct sign */
+        mod = (b > 0) ? 0.0@c@ : -0.0@c@;
+    }
+
+    /* snap quotient to nearest integral value */
+    if (div) {
+        floordiv = npy_floor(div);
+        if (div - floordiv > 0.5@c@)
+            floordiv += 1.0@c@;
+    }
+    else {
+        /* if div is zero ensure correct sign */
+        floordiv = (a / b > 0) ?  0.0@c@ : -0.0@c@;
+    }
+
+    *modulus = mod;
+    return floordiv;
+}
+
 #undef LOGE2
 #undef LOG2E
 #undef RAD2DEG