Fix some things that I was mistakenly testing with an old compiler.

sqrt should not trap on qNaN, which requires some care with
comparisons.  Further, since sqrt(-0) is supposed to be -0, we can't
just use flo:safe-negative? (which returns true for -0.); we must
instead use (flo:safe< x 0.) (which returns false for -0.).

diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm
index ae273999f6cda61d8f06c9f05850e5116cc20df5..7243e2e86e365cb3e1ba743655f6886b093086bc 100644 (file)
--- a/src/runtime/arith.scm
+++ b/src/runtime/arith.scm
@@ -1087,6 +1087,9 @@ USA.
 (define (real:zero? x)
   (if (flonum? x) (flo:zero? x) ((copy rat:zero?) x)))
 
+(define (real:safe-zero? x)
+  (if (flonum? x) (flo:safe-zero? x) ((copy rat:zero?) x)))
+
 (define (real:exact0= x)
   (if (flonum? x) #f ((copy rat:zero?) x)))
 
@@ -1994,19 +1997,19 @@ USA.
   (cond ((recnum? z)
         (let ((x (rec:real-part z))
               (y (rec:imag-part z)))
-          (cond ((real:zero? x)
+          (cond ((real:safe-zero? y)
+                 (assert (not (real:exact0= y)))
+                 (if (if (flonum? x) (flo:safe< x 0.) (rat:negative? x))
+                     (complex:%make-rectangular
+                      0.
+                      (real:copysign (x>=0 (real:negate x)) y))
+                     (complex:%make-rectangular (x>=0 (real:->inexact x)) y)))
+                ((real:safe-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)
-                 (assert (not (real:exact0= y)))
-                 (if (real:negative? x)
-                     (complex:%make-rectangular
-                      0.
-                      (real:copysign (x>=0 (real:negate x)) y))
-                     (complex:%make-rectangular (x>=0 x) y)))
                 ((eq? (real:infinite? x) (real:infinite? y))
                  ;; Standard formula.  Works when both inputs are
                  ;; finite, when both inputs are infinite, and when
@@ -2027,7 +2030,7 @@ USA.
                  ;; NaNity as possible.
                  (assert (or (real:nan? x) (real:nan? y)))
                  (complex:%make-rectangular (flo:abs x) y)))))
-       ((real:negative? z)
+       ((if (flonum? z) (flo:safe< z 0.) (rat:negative? z))
         (complex:%make-rectangular 0 (x>=0 (real:negate z))))
        (else
         (x>=0 z))))

diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm
index 0de60324a8daa3cd514b09d4f66b592ceb204024..82cf23ebfad57e673c6a3a1c7e728dbddad27768 100644 (file)
--- a/tests/runtime/test-arith.scm
+++ b/tests/runtime/test-arith.scm
@@ -748,11 +748,15 @@ USA.
     (,(make-rectangular 0 (- (* 2 (expt 2 4000))))
      ,(make-rectangular (expt 2 2000) (- (expt 2 2000)))
      ,expect-error)
-    ;; Handle signed zero carefully.
+    ;; Handle signed zero carefully.  IEEE 754-2008 specifies that
+    ;; sqrt(-0) = -0, so I guess we'll keep that for the complex
+    ;; extension, but I'm not attached to that.
     (+0.i 0.+0.i)
     (-0.i 0.-0.i)
-    (-0.+0.i +0.+0.i)
-    (-0.-0.i +0.-0.i)
+    (+0.+0.i +0.+0.i)
+    (+0.-0.i +0.-0.i)
+    (-0.+0.i -0.+0.i)
+    (-0.-0.i -0.-0.i)
     ;; Treat infinities carefully around branch cuts.
     (-inf.0 +inf.0i)
     (+inf.0 +inf.0)