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 numpy#12841.

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
 =============

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;

numpy/core/tests/test_umath.py

@@ -273,6 +273,12 @@ def test_floor_division_complex(self):
         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))