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