BUG: Use -0. as initial value for summation

Technically, we should ensure that we do all summations starting with
-0 unless there is nothing to sum (in which case the result is clearly
0).
This is a start, since the initial value is still filled in as 0 by
the reduce machinery.  However, it fixes the odd case where an
inplace addition:

   x1 = np.array(-0.0)
   x2 = np.array(-0.0)
   x1 += x2

May go into the reduce code path (becaus strides are 0).  We could
avoid the reduce path there, but -0 here is strictly correct anyway
and also a necessary step towards fixing `sum([-0., -0., -0.])`
which still requires `initial=-0.0` to be passed manually right now.

There are new `xfail` marked tests also checking the path without
initial.  Presumably, we should be using -0.0 as initial value,
but 0.0 as default (if an array is empty) here.
This may be more reasonably possible after gh-20970.

Closes gh-21213, gh-21212

2 files changed, 61 insertions, 3 deletions

diff --git a/numpy/core/src/umath/loops_utils.h.src b/numpy/core/src/umath/loops_utils.h.src
index 762e9ee59..df92bc315 100644
--- a/numpy/core/src/umath/loops_utils.h.src
+++ b/numpy/core/src/umath/loops_utils.h.src
@@ -79,7 +79,11 @@ static NPY_INLINE @type@
 {
     if (n < 8) {
         npy_intp i;
-        @type@ res = 0.;
+        /*
+         * Start with -0 to preserve -0 values.  The reason is that summing
+         * only -0 should return -0, but `0 + -0 == 0` while `-0 + -0 == -0`.
+         */
+        @type@ res = -0.0;
 
         for (i = 0; i < n; i++) {
             res += @trf@(*((@dtype@*)(a + i * stride)));
@@ -156,8 +160,8 @@ static NPY_INLINE void
     if (n < 8) {
         npy_intp i;
 
-        *rr = 0.;
-        *ri = 0.;
+        *rr = -0.0;
+        *ri = -0.0;
         for (i = 0; i < n; i += 2) {
             *rr += *((@ftype@ *)(a + i * stride + 0));
             *ri += *((@ftype@ *)(a + i * stride + sizeof(@ftype@)));
diff --git a/numpy/core/tests/test_ufunc.py b/numpy/core/tests/test_ufunc.py
index f13c2667d..cb0f7bcbd 100644
--- a/numpy/core/tests/test_ufunc.py
+++ b/numpy/core/tests/test_ufunc.py
@@ -2507,3 +2507,57 @@ def test_ufunc_methods_floaterrors(method):
     with np.errstate(all="raise"):
         with pytest.raises(FloatingPointError):
             method(arr)
+
+
+def _check_neg_zero(value):
+    if value != 0.0:
+        return False
+    if not np.signbit(value.real):
+        return False
+    if value.dtype.kind == "c":
+        return np.signbit(value.imag)
+    return True
+
+@pytest.mark.parametrize("dtype", np.typecodes["AllFloat"])
+def test_addition_negative_zero(dtype):
+    dtype = np.dtype(dtype)
+    if dtype.kind == "c":
+        neg_zero = dtype.type(complex(-0.0, -0.0))
+    else:
+        neg_zero = dtype.type(-0.0)
+
+    arr = np.array(neg_zero)
+    arr2 = np.array(neg_zero)
+
+    assert _check_neg_zero(arr + arr2)
+    # In-place ops may end up on a different path (reduce path) see gh-21211
+    arr += arr2
+    assert _check_neg_zero(arr)
+
+
+@pytest.mark.parametrize("dtype", np.typecodes["AllFloat"])
+@pytest.mark.parametrize("use_initial", [True, False])
+def test_addition_reduce_negative_zero(dtype, use_initial):
+    dtype = np.dtype(dtype)
+    if dtype.kind == "c":
+        neg_zero = dtype.type(complex(-0.0, -0.0))
+    else:
+        neg_zero = dtype.type(-0.0)
+
+    kwargs = {}
+    if use_initial:
+        kwargs["initial"] = neg_zero
+    else:
+        pytest.xfail("-0. propagation in sum currently requires initial")
+
+    # Test various length, in case SIMD paths or chunking play a role.
+    # 150 extends beyond the pairwise blocksize; probably not important.
+    for i in range(0, 150):
+        arr = np.array([neg_zero] * i, dtype=dtype)
+        res = np.sum(arr, **kwargs)
+        if i > 0 or use_initial:
+            assert _check_neg_zero(res)
+        else:
+            # `sum([])` should probably be 0.0 and not -0.0 like `sum([-0.0])`
+            assert not np.signbit(res.real)
+            assert not np.signbit(res.imag)