numpy · charris · Feb 11, 2016 · Feb 12, 2016
diff --git 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_remainder(npy_half x, npy_half y);
 
 /*
  * Half-precision constants

diff --git a/numpy/core/include/numpy/npy_math.h b/numpy/core/include/numpy/npy_math.h
@@ -309,16 +309,19 @@ 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_remainder(double x, double y);
 
 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_remainderf(float x, float y);
 
 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_remainderl(npy_longdouble x, npy_longdouble y);
 
 #define npy_degrees npy_rad2deg
 #define npy_degreesf npy_rad2degf

diff --git a/numpy/core/setup.py b/numpy/core/setup.py
@@ -896,6 +896,7 @@ def generate_umath_c(ext, build_dir):
 
     umath_deps = [
             generate_umath_py,
+            join('include', 'numpy', 'npy_math.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/npymath/halffloat.c b/numpy/core/src/npymath/halffloat.c
@@ -72,6 +72,25 @@ int npy_half_signbit(npy_half h)
     return (h&0x8000u) != 0;
 }
 
+npy_half npy_half_remainder(npy_half x, npy_half y)
+{
+    const npy_half half_zero = (npy_half)0;
+    const float xf = npy_half_to_float(x);
+    const float yf = npy_half_to_float(y);
+    float remf;
+    npy_half remh;
+
+    remh = npy_float_to_half(npy_remainderf(xf, yf));
+    remf = npy_half_to_float(remh);
+    if (yf > 0 && remf >= yf) {
+        remh = npy_half_nextafter(remh, half_zero);
+    }
+    if (yf < 0 && remf <= yf) {
+        remh = npy_half_nextafter(remh, half_zero);
+    }
+    return remh;
+}
+
 npy_half npy_half_spacing(npy_half h)
 {
     npy_half ret;

diff --git a/numpy/core/src/npymath/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src
@@ -608,6 +608,36 @@ double npy_log2(double x)
     }
 }
 
+
+/* remainder(x, y)
+ *
+ * Unlike Python, we assume that the floor function is sacred rather
+ * than fmod. The result is guaranteed to have the same sign as the
+ * divisor and abs(remainder) < abs(y).
+ */
+@type@ npy_remainder@c@(@type@ x, @type@ y)
+{
+    @type@ rem = x - y*npy_floor@c@(x/y);
+
+    if (y < 0) {
+        if (rem >= 0) {
+            rem = -0.0@c@;
+        }
+        else if (rem <= y) {
+            rem = npy_nextafter@c@(y, 0);
+        }
+    }
+    else if (y > 0) {
+        if (rem <= 0) {
+            rem = 0.0@c@;
+        }
+        else if (rem >= y) {
+            rem = npy_nextafter@c@(y, 0);
+        }
+    }
+    return rem;
+}
+
 #undef LOGE2
 #undef LOG2E
 #undef RAD2DEG

diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
@@ -1706,8 +1706,26 @@ 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));
+        @type@ rem;
+
+        rem = in1 - in2*npy_floor@c@(in1/in2);
+        if (in2 < 0) {
+            if (rem >= 0) {
+                rem = -0.0@c@;
+            }
+            else if (rem <= in2) {
+                rem = npy_nextafter@c@(in2, 0);
+            }
+        }
+        else if (in2 > 0) {
+            if (rem <= 0) {
+                rem = 0.0@c@;
+            }
+            else if (rem >= in2) {
+                rem = npy_nextafter@c@(in2, 0);
+            }
+        }
+        *((@type@ *)op1) = rem;
     }
 }
 
@@ -2023,13 +2041,18 @@ HALF_remainder(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNU
     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);
+        float remf;
+        npy_half remh;
+
+        remh = npy_float_to_half(npy_remainderf(in1, in2));
+        remf = npy_half_to_float(remh);
+        if (in2 > 0 && remf >= in2) {
+            remh = npy_half_nextafter(remh, NPY_HALF_ZERO);
         }
-        else {
-            *((npy_half *)op1) = npy_float_to_half(res);
+        if (in2 < 0 && remf <= in2) {
+            remh = npy_half_nextafter(remh, NPY_HALF_ZERO);
         }
+        *((npy_half *)op1) = remh;
     }
 }
 

diff --git a/numpy/core/src/umath/scalarmath.c.src b/numpy/core/src/umath/scalarmath.c.src
@@ -21,6 +21,7 @@
 #include "npy_pycompat.h"
 
 #include "numpy/halffloat.h"
+#include "numpy/npy_math.h"
 #include "templ_common.h"
 
 /* Basic operations:
@@ -283,7 +284,6 @@ static @type@ (*_basic_@name@_fmod)(@type@, @type@);
 
 static npy_half (*_basic_half_floor)(npy_half);
 static npy_half (*_basic_half_sqrt)(npy_half);
-static npy_half (*_basic_half_fmod)(npy_half, npy_half);
 
 #define half_ctype_add(a, b, outp) *(outp) = \
         npy_float_to_half(npy_half_to_float(a) + npy_half_to_float(b))
@@ -353,22 +353,21 @@ static npy_half (*_basic_half_fmod)(npy_half, npy_half);
     } while(0)
 /**end repeat**/
 
+
 /**begin repeat
  * #name = float, double, longdouble#
  * #type = npy_float, npy_double, npy_longdouble#
+ * #c = f, ,l#
  */
 static void
 @name@_ctype_remainder(@type@ a, @type@ b, @type@ *out) {
-    @type@ tmp = a/b;
-    *out = b * (tmp - _basic_@name@_floor(tmp));
+    *out = npy_remainder@c@(a, b);
 }
 /**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);
+    *out = npy_half_remainder(a, b);
 }
 
 
@@ -1721,7 +1720,6 @@ get_functions(PyObject * mm)
         i += 3;
         j++;
     }
-    _basic_half_fmod = funcdata[j - 1];
     _basic_float_fmod = funcdata[j];
     _basic_double_fmod = funcdata[j + 1];
     _basic_longdouble_fmod = funcdata[j + 2];

diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
@@ -2448,42 +2448,6 @@ def test_conjugate(self):
         assert_raises(AttributeError, lambda: a.conj())
         assert_raises(AttributeError, lambda: a.conjugate())
 
-    def test_divmod_basic(self):
-        dt = np.typecodes['AllInteger'] + np.typecodes['Float']
-        for dt1, dt2 in itertools.product(dt, dt):
-            for sg1, sg2 in itertools.product((+1, -1), (+1, -1)):
-                if sg1 == -1 and dt1 in np.typecodes['UnsignedInteger']:
-                    continue
-                if sg2 == -1 and dt2 in np.typecodes['UnsignedInteger']:
-                    continue
-                fmt = 'dt1: %s, dt2: %s, sg1: %s, sg2: %s'
-                msg = fmt % (dt1, dt2, sg1, sg2)
-                a = np.array(sg1*71, dtype=dt1)
-                b = np.array(sg2*19, dtype=dt2)
-                div, rem = divmod(a, b)
-                assert_allclose(div*b + rem, a, err_msg=msg)
-                if sg2 == -1:
-                    assert_(b < rem <= 0, msg)
-                else:
-                    assert_(b > rem >= 0, msg)
-
-    def test_divmod_roundoff(self):
-        # gh-6127
-        dt = 'fdg'
-        for dt1, dt2 in itertools.product(dt, dt):
-            for sg1, sg2 in itertools.product((+1, -1), (+1, -1)):
-                fmt = 'dt1: %s, dt2: %s, sg1: %s, sg2: %s'
-                msg = fmt % (dt1, dt2, sg1, sg2)
-                a = np.array(sg1*78*6e-8, dtype=dt1)
-                b = np.array(sg2*6e-8, dtype=dt2)
-                div, rem = divmod(a, b)
-                assert_allclose(div*b + rem, a, err_msg=msg)
-                if sg2 == -1:
-                    assert_(b < rem <= 0, msg)
-                else:
-                    assert_(b > rem >= 0, msg)
-
-
 class TestBinop(object):
     def test_inplace(self):
         # test refcount 1 inplace conversion

diff --git a/numpy/core/tests/test_scalarmath.py b/numpy/core/tests/test_scalarmath.py
@@ -2,6 +2,8 @@
 
 import sys
 import itertools
+import warnings
+import operator
 
 import numpy as np
 from numpy.testing.utils import _gen_alignment_data
@@ -137,8 +139,12 @@ def test_mixed_types(self):
                     assert_almost_equal(result, 9, err_msg=msg)
 
 
-class TestDivmod(TestCase):
-    def test_divmod_basic(self):
+class TestModulus(TestCase):
+
+    floordiv = operator.floordiv
+    mod = operator.mod
+
+    def test_modulus_basic(self):
         dt = np.typecodes['AllInteger'] + np.typecodes['Float']
         for dt1, dt2 in itertools.product(dt, dt):
             for sg1, sg2 in itertools.product((+1, -1), (+1, -1)):
@@ -150,29 +156,86 @@ def test_divmod_basic(self):
                 msg = fmt % (dt1, dt2, sg1, sg2)
                 a = np.array(sg1*71, dtype=dt1)[()]
                 b = np.array(sg2*19, dtype=dt2)[()]
-                div, rem = divmod(a, b)
+                div = self.floordiv(a, b)
+                rem = self.mod(a, b)
                 assert_allclose(div*b + rem, a, err_msg=msg)
                 if sg2 == -1:
                     assert_(b < rem <= 0, msg)
                 else:
                     assert_(b > rem >= 0, msg)
 
-    def test_divmod_roundoff(self):
+    def test_float_modulus_exact(self):
+        # test that float results are exact for small integers. This also
+        # holds for the same integers scaled by powers of two.
+        nlst = list(range(-127, 0))
+        plst = list(range(1, 128))
+        dividend = nlst + [0] + plst
+        divisor = nlst + plst
+        arg = list(itertools.product(dividend, divisor))
+        tgt = list(divmod(*t) for t in arg)
+
+        a, b = np.array(arg, dtype=int).T
+        # convert exact integer results from Python to float so that
+        # signed zero can be used, it is checked.
+        tgtdiv, tgtrem = np.array(tgt, dtype=float).T
+        tgtdiv = np.where((tgtdiv == 0.0) & ((b < 0) ^ (a < 0)), -0.0, tgtdiv)
+        tgtrem = np.where((tgtrem == 0.0) & (b < 0), -0.0, tgtrem)
+
+        for dt in np.typecodes['Float']:
+            msg = 'dtype: %s' % (dt,)
+            fa = a.astype(dt)
+            fb = b.astype(dt)
+            # use list comprehension so a_ and b_ are scalars
+            div = [self.floordiv(a_, b_) for  a_, b_ in zip(fa, fb)]
+            rem = [self.mod(a_, b_) for a_, b_ in zip(fa, fb)]
+            assert_equal(div, tgtdiv, err_msg=msg)
+            assert_equal(rem, tgtrem, err_msg=msg)
+
+    def test_float_modulus_roundoff(self):
         # gh-6127
-        dt = 'fdg'
+        dt = np.typecodes['Float']
         for dt1, dt2 in itertools.product(dt, dt):
             for sg1, sg2 in itertools.product((+1, -1), (+1, -1)):
                 fmt = 'dt1: %s, dt2: %s, sg1: %s, sg2: %s'
                 msg = fmt % (dt1, dt2, sg1, sg2)
                 a = np.array(sg1*78*6e-8, dtype=dt1)[()]
                 b = np.array(sg2*6e-8, dtype=dt2)[()]
-                div, rem = divmod(a, b)
+                div = self.floordiv(a, b)
+                rem = self.mod(a, b)
                 assert_allclose(div*b + rem, a, err_msg=msg)
                 if sg2 == -1:
                     assert_(b < rem <= 0, msg)
                 else:
                     assert_(b > rem >= 0, msg)
 
+    def test_float_modulus_corner_cases(self):
+        # Check remainder magnitude.
+        for dt in np.typecodes['Float']:
+            b = np.array(1.0, dtype=dt)
+            a = np.nextafter(np.array(0.0, dtype=dt), -b)
+            rem = self.mod(a, b)
+            assert_(rem < b, 'dt: %s' % dt)
+            rem = self.mod(-a, -b)
+            assert_(rem > -b, 'dt: %s' % dt)
+
+        # Check nans, inf
+        with warnings.catch_warnings():
+            warnings.simplefilter('always')
+            warnings.simplefilter('ignore', RuntimeWarning)
+            for dt in np.typecodes['Float']:
+                fone = np.array(1.0, dtype=dt)
+                fzer = np.array(0.0, dtype=dt)
+                finf = np.array(np.inf, dtype=dt)
+                fnan = np.array(np.nan, dtype=dt)
+                rem = self.mod(fone, fzer)
+                assert_(np.isnan(rem), 'dt: %s' % dt)
+                rem = self.mod(fone, finf)
+                assert_(np.isnan(rem), 'dt: %s' % dt)
+                rem = self.mod(fone, fnan)
+                assert_(np.isnan(rem), 'dt: %s' % dt)
+                rem = self.mod(finf, fone)
+                assert_(np.isnan(rem), 'dt: %s' % dt)
+
 
 class TestComplexDivision(TestCase):
     def test_zero_division(self):