[src/trunk]: src/lib/libm/complex csqrt has a branch cut on the negative real...

[maya <maya%NetBSD.org@localhost>]

details:   https://anonhg.NetBSD.org/src/rev/7f49c0b9f3ae
branches:  trunk
changeset: 349998:7f49c0b9f3ae
user:      maya <maya%NetBSD.org@localhost>
date:      Sat Dec 31 20:01:15 2016 +0000

description:
csqrt has a branch cut on the negative real axis, and this requires
delicacy in order to maintain continuity around it.

we have an initial case to deal with a fairly common case: getting
a real number. Avoid dealing with the branch cut in this case by
checking if the real part is negative.

later, -0.0 < 0 is not met, so instead, test for a negative number
using signbit, so negative zero is also treated as a negative number.

Fixes last part of PR lib/51427: libm issues triggered by py-numpy

ok riastradh

diffstat:

 lib/libm/complex/csqrt.c  |  7 ++++---
 lib/libm/complex/csqrtf.c |  7 ++++---
 2 files changed, 8 insertions(+), 6 deletions(-)

diffs (56 lines):

diff -r f1ef64512ba8 -r 7f49c0b9f3ae lib/libm/complex/csqrt.c
--- a/lib/libm/complex/csqrt.c  Sat Dec 31 17:46:35 2016 +0000
+++ b/lib/libm/complex/csqrt.c  Sat Dec 31 20:01:15 2016 +0000
@@ -1,4 +1,4 @@
-/* $NetBSD: csqrt.c,v 1.2 2016/12/31 15:33:03 maya Exp $ */
+/* $NetBSD: csqrt.c,v 1.3 2016/12/31 20:01:15 maya Exp $ */
 
 /*-
  * Copyright (c) 2007 The NetBSD Foundation, Inc.
@@ -41,7 +41,8 @@
        x = creal (z);
        y = cimag (z);
 
-       if (y == 0.0) {
+       /* Input is a real number that isn't on the branch cut */
+       if ((y == 0.0) && !signbit(y)) {
                if (x == 0.0) {
                        w = 0.0 + y * I;
                } else {
@@ -92,7 +93,7 @@
                t = scale * fabs((0.5 * y) / r);
                r *= scale;
        }
-       if (y < 0)
+       if (signbit(y))
                w = t - r * I;
        else
                w = t + r * I;
diff -r f1ef64512ba8 -r 7f49c0b9f3ae lib/libm/complex/csqrtf.c
--- a/lib/libm/complex/csqrtf.c Sat Dec 31 17:46:35 2016 +0000
+++ b/lib/libm/complex/csqrtf.c Sat Dec 31 20:01:15 2016 +0000
@@ -1,4 +1,4 @@
-/* $NetBSD: csqrtf.c,v 1.1 2007/08/20 16:01:37 drochner Exp $ */
+/* $NetBSD: csqrtf.c,v 1.2 2016/12/31 20:01:15 maya Exp $ */
 
 /*-
  * Copyright (c) 2007 The NetBSD Foundation, Inc.
@@ -41,7 +41,8 @@
        x = crealf (z);
        y = cimagf (z);
 
-       if (y == 0.0f) {
+       /* Input is a real number that isn't on the branch cut */
+       if ((y == 0.0f) && !signbit(y)) {
                if (x < 0.0f) {
                        w = 0.0f + sqrtf(-x) * I;
                        return w;
@@ -91,7 +92,7 @@
                r *= scale;
        }
 
-       if (y < 0)
+       if (signbit(y))
                w = t - r * I;
        else
                w = t + r * I;