BUG: fix signed zero behavior in npy_divmod

Previously when doing floor division numpy would sometimes return an incorrect
signed zero. For example:

>>> np.zeros(10)//1
array([-0., -0., -0., -0., -0., -0., -0., -0., -0., -0.])

>>> np.remainder(-1.0,1.0)
-0.0

The reason for this is that whenever div or mod were zero the code was using
the following to pick the sign of zero:

floordiv = (a / b > 0) ? 0.0@c@ : -0.0@c@;

This commit updates these lines to instead use the copysign function which is
how cpython does floor division.

Fixes #12841.

3 files changed, 22 insertions, 2 deletions

diff --git a/doc/release/1.17.0-notes.rst b/doc/release/1.17.0-notes.rst
index d3a7779c8..57f7505ff 100644
--- a/doc/release/1.17.0-notes.rst
+++ b/doc/release/1.17.0-notes.rst
@@ -36,7 +36,21 @@ Casting from a different floating point precision to float16 used incorrect
 rounding in some edge cases. This means in rare cases, subnormal results will
 now be rounded up instead of down, changing the last bit (ULP) of the result.
 
+Signed zero when using divmod
+-----------------------------
 
+Starting in version 1.12.0, numpy incorrectly returned a negatively signed zero
+when using the ``divmod`` and ``floor_divide`` functions when the result was
+zero. For example::
+
+   >>> np.zeros(10)//1
+   array([-0., -0., -0., -0., -0., -0., -0., -0., -0., -0.])
+
+With this release, the result is correctly returned as a positively signed
+zero::
+
+   >>> np.zeros(10)//1
+   array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
 
 C API changes
 =============
diff --git a/numpy/core/src/npymath/npy_math_internal.h.src b/numpy/core/src/npymath/npy_math_internal.h.src
index f2e5229b0..fa820baac 100644
--- a/numpy/core/src/npymath/npy_math_internal.h.src
+++ b/numpy/core/src/npymath/npy_math_internal.h.src
@@ -654,7 +654,7 @@ npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus)
     }
     else {
         /* if mod is zero ensure correct sign */
-        mod = (b > 0) ? 0.0@c@ : -0.0@c@;
+        mod = npy_copysign@c@(0, b);
     }
 
     /* snap quotient to nearest integral value */
@@ -665,7 +665,7 @@ npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus)
     }
     else {
         /* if div is zero ensure correct sign */
-        floordiv = (a / b > 0) ?  0.0@c@ : -0.0@c@;
+        floordiv = npy_copysign@c@(0, a/b);
     }
 
     *modulus = mod;
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 21097244f..db6b51922 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -273,6 +273,12 @@ class TestDivision(object):
         y = np.floor_divide(x**2, x)
         assert_equal(y, [1.e+110, 0], err_msg=msg)
 
+    def test_floor_division_signed_zero(self):
+        # Check that the sign bit is correctly set when dividing positive and
+        # negative zero by one.
+        x = np.zeros(10)
+        assert_equal(np.signbit(x//1), 0)
+        assert_equal(np.signbit((-x)//1), 1)
 
 def floor_divide_and_remainder(x, y):
     return (np.floor_divide(x, y), np.remainder(x, y))