ENH: Make numpy floating scalars consistent with Python divmod.

The following numpy scalar floating functions are reimplemented
using the npy_divmod function.

- remainder ('%')
- floor_division ('//')
- divmod

Note that numpy scalars warn on zero division rather than raise.

divmod example, Python floats

    In [1]: a = 78 * 6e-8

    In [2]: b = 6e-8

    In [3]: divmod(a, b)
    Out[3]: (77.0, 5.999999999999965e-08)

Before this commit numpy float64 gave

    In [4]: divmod(float64(a), float64(b))
    Out[4]: (78.0, 0.0)

After this commit numpy float64 gives

    In [4]: divmod(float64(a), float64(b))
    Out[4]: (77.0, 5.9999999999999651e-08)

2 files changed, 55 insertions, 27 deletions

diff --git a/numpy/core/setup.py b/numpy/core/setup.py
index 2e9e277af..593b21f6d 100644
--- a/numpy/core/setup.py
+++ b/numpy/core/setup.py
@@ -896,6 +896,8 @@ def configuration(parent_package='',top_path=None):
 
     umath_deps = [
             generate_umath_py,
+            join('include', 'numpy', 'npy_math.h'),
+            join('include', 'numpy', 'halffloat.h'),
             join('src', 'multiarray', 'common.h'),
             join('src', 'private', 'templ_common.h.src'),
             join('src', 'umath', 'simd.inc.src'),
diff --git a/numpy/core/src/umath/scalarmath.c.src b/numpy/core/src/umath/scalarmath.c.src
index 706eccb31..0f762d707 100644
--- a/numpy/core/src/umath/scalarmath.c.src
+++ b/numpy/core/src/umath/scalarmath.c.src
@@ -175,6 +175,7 @@ static void
 }
 
 #define @name@_ctype_floor_divide @name@_ctype_divide
+
 static void
 @name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) {
     if (a == 0 || b == 0) {
@@ -268,17 +269,39 @@ static void
 /**begin repeat
  * #name = float, double, longdouble#
  * #type = npy_float, npy_double, npy_longdouble#
+ * #c = f, , l#
  */
 static @type@ (*_basic_@name@_floor)(@type@);
 static @type@ (*_basic_@name@_sqrt)(@type@);
 static @type@ (*_basic_@name@_fmod)(@type@, @type@);
+
 #define @name@_ctype_add(a, b, outp) *(outp) = (a) + (b)
 #define @name@_ctype_subtract(a, b, outp) *(outp) = (a) - (b)
 #define @name@_ctype_multiply(a, b, outp) *(outp) = (a) * (b)
 #define @name@_ctype_divide(a, b, outp) *(outp) = (a) / (b)
 #define @name@_ctype_true_divide @name@_ctype_divide
-#define @name@_ctype_floor_divide(a, b, outp)   \
-    *(outp) = _basic_@name@_floor((a) / (b))
+
+
+static void
+@name@_ctype_floor_divide(@type@ a, @type@ b, @type@ *out) {
+    @type@ mod;
+
+    *out = npy_divmod@c@(a, b, &mod);
+}
+
+
+static void
+@name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) {
+    npy_divmod@c@(a, b, out);
+}
+
+
+static void
+@name@_ctype_divmod(@type@ a, @type@ b, @type@ *out1, @type@ *out2) {
+    *out1 = npy_divmod@c@(a, b, out2);
+}
+
+
 /**end repeat**/
 
 static npy_half (*_basic_half_floor)(npy_half);
@@ -294,9 +317,26 @@ static npy_half (*_basic_half_fmod)(npy_half, npy_half);
 #define half_ctype_divide(a, b, outp) *(outp) = \
         npy_float_to_half(npy_half_to_float(a) / npy_half_to_float(b))
 #define half_ctype_true_divide half_ctype_divide
-#define half_ctype_floor_divide(a, b, outp)   \
-        *(outp) = npy_float_to_half(_basic_float_floor(npy_half_to_float(a) / \
-        npy_half_to_float(b)))
+
+
+static void
+half_ctype_floor_divide(npy_half a, npy_half b, npy_half *out) {
+    npy_half mod;
+
+    *out = npy_half_divmod(a, b, &mod);
+}
+
+
+static void
+half_ctype_remainder(npy_half a, npy_half b, npy_half *out) {
+    npy_half_divmod(a, b, out);
+}
+
+
+static void
+half_ctype_divmod(npy_half a, npy_half b, npy_half *out1, npy_half *out2) {
+    *out1 = npy_half_divmod(a, b, out2);
+}
 
 /**begin repeat
  * #name = cfloat, cdouble, clongdouble#
@@ -344,38 +384,23 @@ static npy_half (*_basic_half_fmod)(npy_half, npy_half);
         (outp)->imag = (in1i*rat - in1r)*scl;                       \
     }                                                               \
     } while(0)
+
 #define @name@_ctype_true_divide @name@_ctype_divide
+
 #define @name@_ctype_floor_divide(a, b, outp) do {      \
-    (outp)->real = _basic_@rname@_floor                 \
-    (((a).real*(b).real + (a).imag*(b).imag) /          \
-     ((b).real*(b).real + (b).imag*(b).imag));          \
+    @rname@_ctype_floor_divide(                         \
+        ((a).real*(b).real + (a).imag*(b).imag),        \
+        ((b).real*(b).real + (b).imag*(b).imag),        \
+        &((outp)->real));                               \
     (outp)->imag = 0;                                   \
     } while(0)
 /**end repeat**/
 
-/**begin repeat
- * #name = float, double, longdouble#
- * #type = npy_float, npy_double, npy_longdouble#
- */
-static void
-@name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) {
-    @type@ tmp = a/b;
-    *out = b * (tmp - _basic_@name@_floor(tmp));
-}
-/**end repeat**/
-
-static void
-half_ctype_remainder(npy_half a, npy_half b, npy_half *out) {
-    float tmp, fa = npy_half_to_float(a), fb = npy_half_to_float(b);
-    float_ctype_remainder(fa, fb, &tmp);
-    *out = npy_float_to_half(tmp);
-}
 
 
 /**begin repeat
  * #name = byte, ubyte, short, ushort, int, uint, long, ulong,
- *         longlong, ulonglong, half, float, double, longdouble,
- *         cfloat, cdouble, clongdouble#
+ *         longlong, ulonglong, cfloat, cdouble, clongdouble#
  */
 #define @name@_ctype_divmod(a, b, out, out2) {  \
     @name@_ctype_floor_divide(a, b, out);       \
@@ -383,6 +408,7 @@ half_ctype_remainder(npy_half a, npy_half b, npy_half *out) {
     }
 /**end repeat**/
 
+
 /**begin repeat
  * #name = float, double, longdouble#
  * #type = npy_float, npy_double, npy_longdouble#