Fix sqrt along negative real axis with inexact zero imaginary part.

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index c620b906838fc322cdf38cc51dce688754a9c7f8..0a50128cc9037daa26c3f4074ff45f777c0f96a0 100644 (file)
--- 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 b4c0191b53e4f2c6cd35221321c27ff437981324..0c449e91c016dfae44aa5051cca54481e758298b 100644 (file)
--- 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 ()