diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py
index e3600406c117..a017ca26d43e 100644
--- a/numpy/core/code_generators/ufunc_docstrings.py
+++ b/numpy/core/code_generators/ufunc_docstrings.py
@@ -1225,8 +1225,10 @@ def add_newdoc(place, name, doc):
 
 add_newdoc('numpy.core.umath', 'floor_divide',
     """
-    Return the largest integer smaller or equal to the division of the
-    inputs.
+    Return the largest integer smaller or equal to the division of the inputs.
+    It is equivalent to the Python ``//`` operator and pairs with the
+    Python ``%`` (`remainder`), function so that ``b = a % b + b * (a // b)``
+    up to roundoff.
 
     Parameters
     ----------
@@ -1243,6 +1245,7 @@ def add_newdoc(place, name, doc):
 
     See Also
     --------
+    remainder : Remainder complementary to floor_divide.
     divide : Standard division.
     floor : Round a number to the nearest integer toward minus infinity.
     ceil : Round a number to the nearest integer toward infinity.
@@ -2689,9 +2692,9 @@ def add_newdoc(place, name, doc):
     """
     Return element-wise remainder of division.
 
-    Computes ``x1 - floor(x1 / x2) * x2``, the result has the same sign as
-    the divisor `x2`. It is equivalent to the Python modulus operator
-    ``x1 % x2`` and should not be confused with the Matlab(TM) ``rem``
+    Computes the remainder complementary to the `floor_divide` function.  It is
+    equivalent to the Python modulus operator``x1 % x2`` and has the same sign
+    as the divisor `x2`. It should not be confused with the Matlab(TM) ``rem``
     function.
 
     Parameters
@@ -2707,11 +2710,12 @@ def add_newdoc(place, name, doc):
     Returns
     -------
     y : ndarray
-        The remainder of the quotient ``x1/x2``, element-wise. Returns a
-        scalar if both  `x1` and `x2` are scalars.
+        The element-wise remainder of the quotient ``floor_divide(x1, x2)``.
+        Returns a scalar if both  `x1` and `x2` are scalars.
 
     See Also
     --------
+    floor_divide : Equivalent of Python ``//`` operator.
     fmod : Equivalent of the Matlab(TM) ``rem`` function.
     divide, floor
 
diff --git a/numpy/core/include/numpy/halffloat.h b/numpy/core/include/numpy/halffloat.h
index 944f0ea34b48..ab0d221fb431 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 3dae583f3f8f..e76508de0479 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/setup.py b/numpy/core/setup.py
index 2e9e277af01c..593b21f6d864 100644
--- a/numpy/core/setup.py
+++ b/numpy/core/setup.py
@@ -896,6 +896,8 @@ def generate_umath_c(ext, build_dir):
 
     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/npymath/halffloat.c b/numpy/core/src/npymath/halffloat.c
index 34ac642876eb..951768256623 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 4dcb01986d7b..45b618a566f7 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
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 7b8dcdbafd6e..0d9806f5d6ab 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);
     }
 }
 
diff --git a/numpy/core/src/umath/scalarmath.c.src b/numpy/core/src/umath/scalarmath.c.src
index 706eccb31761..77520abf557f 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,20 +269,40 @@ 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);
 static npy_half (*_basic_half_sqrt)(npy_half);
 static npy_half (*_basic_half_fmod)(npy_half, npy_half);
 
@@ -294,9 +315,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 +382,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 +406,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#
@@ -1665,25 +1689,6 @@ get_functions(PyObject * mm)
     _basic_clongdouble_pow = funcdata[j + 5];
     Py_DECREF(obj);
 
-    /* Get the floor functions */
-    obj = PyObject_GetAttrString(mm, "floor");
-    if (obj == NULL) {
-        goto fail;
-    }
-    funcdata = ((PyUFuncObject *)obj)->data;
-    signatures = ((PyUFuncObject *)obj)->types;
-    i = 0;
-    j = 0;
-    while(signatures[i] != NPY_FLOAT) {
-        i += 2;
-        j++;
-    }
-    _basic_half_floor = funcdata[j - 1];
-    _basic_float_floor = funcdata[j];
-    _basic_double_floor = funcdata[j + 1];
-    _basic_longdouble_floor = funcdata[j + 2];
-    Py_DECREF(obj);
-
     /* Get the sqrt functions */
     obj = PyObject_GetAttrString(mm, "sqrt");
     if (obj == NULL) {
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index d57e7c106751..16df605f4fa7 100644
--- 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
index 17f70f6c96a8..12b1a0fe3356 100644
--- 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,88 @@ 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)
-                assert_allclose(div*b + rem, a, err_msg=msg)
+                div = self.floordiv(a, b)
+                rem = self.mod(a, b)
+                assert_equal(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)
-                assert_allclose(div*b + rem, a, err_msg=msg)
+                div = self.floordiv(a, b)
+                rem = self.mod(a, b)
+                # Equal assertion should hold when fmod is used
+                assert_equal(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)
+                # MSVC 2008 returns NaN here, so disable the check.
+                #rem = self.mod(fone, finf)
+                #assert_(rem == fone, '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):
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 917e05e6a3c4..da52e0dde394 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -3,6 +3,7 @@
 import sys
 import platform
 import warnings
+import itertools
 
 from numpy.testing.utils import _gen_alignment_data
 import numpy.core.umath as ncu
@@ -222,6 +223,102 @@ def test_floor_division_complex(self):
         assert_equal(y, [1.e+110, 0], err_msg=msg)
 
 
+class TestRemainder(TestCase):
+
+    def test_remainder_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 = np.floor_divide(a, b)
+                rem = np.remainder(a, b)
+                assert_equal(div*b + rem, a, err_msg=msg)
+                if sg2 == -1:
+                    assert_(b < rem <= 0, msg)
+                else:
+                    assert_(b > rem >= 0, msg)
+
+    def test_float_remainder_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)
+            div = np.floor_divide(fa, fb)
+            rem = np.remainder(fa, fb)
+            assert_equal(div, tgtdiv, err_msg=msg)
+            assert_equal(rem, tgtrem, err_msg=msg)
+
+    def test_float_remainder_roundoff(self):
+        # gh-6127
+        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 = np.floor_divide(a, b)
+                rem = np.remainder(a, b)
+                # Equal assertion should hold when fmod is used
+                assert_equal(div*b + rem, a, err_msg=msg)
+                if sg2 == -1:
+                    assert_(b < rem <= 0, msg)
+                else:
+                    assert_(b > rem >= 0, msg)
+
+    def test_float_remainder_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 = np.remainder(a, b)
+            assert_(rem <= b, 'dt: %s' % dt)
+            rem = np.remainder(-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 = np.remainder(fone, fzer)
+                assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
+                # MSVC 2008 returns NaN here, so disable the check.
+                #rem = np.remainder(fone, finf)
+                #assert_(rem == fone, 'dt: %s, rem: %s' % (dt, rem))
+                rem = np.remainder(fone, fnan)
+                assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
+                rem = np.remainder(finf, fone)
+                assert_(np.isnan(rem), 'dt: %s, rem: %s' % (dt, rem))
+
+
 class TestCbrt(TestCase):
     def test_cbrt_scalar(self):
         assert_almost_equal((np.cbrt(np.float32(-2.5)**3)), -2.5)