Preserve signed-zero branch sides in complex square roots

RELEASE_NOTES.rst

@@ -6,6 +6,8 @@ Changes from 2.14.2 to 2.14.3
 -----------------------------
 
 * **Under development.**
+* Preserve the sign of the imaginary component in complex square roots on
+  the negative real axis, including signed zero inputs.
 * Replaced expression-string blacklist filtering with AST validation before
   evaluation, and disabled Python builtins while sanitization is enabled.
   This closes sanitizer bypasses in both cached and ``disable_cache=True``

numexpr/complex_functions.hpp

@@ -109,22 +109,20 @@ static void
 nc_sqrt(std::complex<double> *x, std::complex<double> *r)
 {
     double s,d;
-    if (x->real() == 0. && x->imag() == 0.)
-        *r = *x;
+    if (x->real() == 0. && x->imag() == 0.) {
+        r->imag(x->imag());
+        r->real(0.);
+    }
     else {
         s = sqrt((fabs(x->real()) + hypot(x->real(),x->imag()))/2);
         d = x->imag()/(2*s);
         if (x->real() > 0.) {
             r->real(s);
             r->imag(d);
         }
-        else if (x->imag() >= 0.) {
-            r->real(d);
-            r->imag(s);
-        }
         else {
-            r->real(-d);
-            r->imag(-s);
+            r->real(fabs(d));
+            r->imag(copysign(s, x->imag()));
         }
     }
     return;

numexpr/tests/test_numexpr.py

@@ -532,6 +532,22 @@ def test_bitwise_operators(self):
         assert_array_equal(evaluate("x | y"), x | y) # or
         assert_array_equal(evaluate("~x"), ~x) # invert
 
+    def test_complex_sqrt_branch_side(self):
+        # Signed imaginary zero distinguishes the two sides of the branch cut.
+        x = np.array([complex(real, imag)
+                      for real in (-4., -0., 0., 4.)
+                      for imag in (-0., 0., -1., 1.)])
+        for values in (x, x[::-1], np.tile(x, 1024)):
+            expected = np.sqrt(values)
+            actual = evaluate("sqrt(values)")
+            assert_allclose(actual, expected)
+            assert_array_equal(np.signbit(actual.real), np.signbit(expected.real))
+            assert_array_equal(np.signbit(actual.imag), np.signbit(expected.imag))
+            assert_allclose(evaluate("sqrt(conj(values))"), np.conj(expected))
+            out = values.copy()
+            evaluate("sqrt(out)", out=out)
+            assert_allclose(out, expected)
+
     def test_complex_tan(self):
         # old version of NumExpr had overflow problems
         x = np.arange(1, 400., step=16., dtype=np.complex128)