From 2a084600e9670127c6b31b32ad71c2789d7c71dd Mon Sep 17 00:00:00 2001
From: Taylor R Campbell <campbell@mumble.net>
Date: Fri, 30 Nov 2018 04:29:55 +0000
Subject: [PATCH] Fix sqrt along negative real axis with inexact zero imaginary
 part.

---
 src/runtime/arith.scm        | 29 +++++++++++++++++++----------
 tests/runtime/test-arith.scm | 16 ++++++++--------
 2 files changed, 27 insertions(+), 18 deletions(-)

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index c620b9068..0a50128cc 100644
--- a/src/runtime/arith.scm
+++ b/src/runtime/arith.scm
@@ -1967,21 +1967,30 @@ USA.
       (real:abs z)))
 
 (define (complex:sqrt z)
+  (define (x>=0 x)
+    ((copy real:sqrt) x))
   (cond ((recnum? z)
 	 (let ((x (rec:real-part z))
 	       (y (rec:imag-part z)))
-	   (if (real:zero? x)
-	       ;; sqrt(+/- 2i x) = sqrt(x) +/- sqrt(x)i
-	       (let ((sqrt-abs-y/2 (real:sqrt (real:/ (real:abs y) 2))))
-		 (complex:%make-rectangular
-		  sqrt-abs-y/2
-		  (real:copysign sqrt-abs-y/2 y)))
-	       (complex:%make-polar (real:sqrt (complex:magnitude z))
-                                    (real:/ (complex:angle z) 2)))))
+	   (cond ((real:zero? x)
+		  ;; sqrt(+/- 2i x) = sqrt(x) +/- i sqrt(x)
+		  (let ((sqrt-abs-y/2 (x>=0 (real:/ (real:abs y) 2))))
+		    (complex:%make-rectangular
+		     sqrt-abs-y/2
+		     (real:copysign sqrt-abs-y/2 y))))
+		 ((real:zero? y)
+		  (if (real:negative? x)
+		      (complex:%make-rectangular
+		       0.
+		       (real:copysign (x>=0 (real:negate x)) y))
+		      (complex:%make-rectangular (x>=0 x) y)))
+		 (else
+		  (complex:%make-polar (x>=0 (complex:magnitude z))
+				       (real:/ (complex:angle z) 2))))))
 	((real:negative? z)
-	 (complex:%make-rectangular 0 (real:sqrt (real:negate z))))
+	 (complex:%make-rectangular 0 (x>=0 (real:negate z))))
 	(else
-	 ((copy real:sqrt) z))))
+	 (x>=0 z))))
 
 (define (complex:expt z1 z2)
   (cond ((complex:zero? z1)
diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index b4c0191b5..0c449e91c 100644
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -714,17 +714,17 @@ USA.
     (-0.i 0.-0.i)
     (-0.+0.i +0.+0.i)
     (-0.-0.i +0.-0.i)
-    (-4.+0.i 0.+2.i ,expect-failure)
-    (-4.-0.i 0.-2.i ,expect-failure)
+    (-4.+0.i 0.+2.i)
+    (-4.-0.i 0.-2.i)
     ;; Treat infinities carefully around branch cuts.
     (-inf.0 +inf.0i)
     (+inf.0 +inf.0)
-    (-inf.0+0.i 0.+inf.0i ,expect-failure)
-    (+inf.0+0.i +inf.0+0.i ,expect-error)
+    (-inf.0+0.i 0.+inf.0i)
+    (+inf.0+0.i +inf.0+0.i)
     (-inf.0+1.i 0.+inf.0i ,expect-failure)
     (+inf.0+1.i +inf.0+0.i ,expect-error)
-    (-inf.0-0.i 0.-inf.0i ,expect-failure)
-    (+inf.0-0.i +inf.0-0.i ,expect-error)
+    (-inf.0-0.i 0.-inf.0i)
+    (+inf.0-0.i +inf.0-0.i)
     (-inf.0-1.i 0.-inf.0i ,expect-failure)
     (+inf.0-1.i +inf.0-0.i ,expect-error)
     (-inf.0i +inf.0-inf.0i)
@@ -747,8 +747,8 @@ USA.
     (-4. +2.i)
     ;; Square root of negative real with inexact zero imaginary part
     ;; should be imaginary with inexact zero real part.
-    (-4.+0.i 0.+2.i ,expect-failure)
-    (-4.-0.i 0.-2.i ,expect-failure))
+    (-4.+0.i 0.+2.i)
+    (-4.-0.i 0.-2.i))
   (lambda (z r #!optional xfail)
     (with-expected-failure xfail
       (lambda ()
-- 
2.25.1