On 2012-12-03, Raymond Toy <toy.raymond at gmail.com> wrote:
> Robert> I dunno -- isn't z = -i a counterexample to that?
>
> See my reply to Jaime. With signed zeroes z = -i = +0 - i, so things
> work out differently and are easier to reason about, at least in this
> case.
Without bringing signed zero into the picture, can you show that
sqrt((-i)^2)/(-i) is something other than -1 ?
Given erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2, z^2)/sqrt(pi)),
if sqrt(z^2)/z is 1 for pure imaginary z, then erf(-z) = erf(z) (which
is how I bumped into this to begin with). I don't think we want that.
> CL doesn't forbid signed zeroes.
Right, but it doesn't require them, so any code depending on them is
nonportable.
For the record, I've committed and pushed 8f10856b70 to address the
problems with erf & friends.
best
Robert Dodier
PS.
commit 8f10856b70ef40660bd75461352ca9ace2d58f2f
Author: robert_dodier <robert_dodier at users.sourceforge.net>
Date: Mon Dec 3 00:26:34 2012 -0700
Work around bugs in COMPLEX-ERF. Bug #1: sqrt(z^2)/z incorrect
when Re(z) = 0 and Im(z) < 0. Bug #2: (EXPT Z 2.0) yields spurious imaginary part
in some Lisps. Bug #3: GAMMA-INCOMPLETE returns different values (conjugates) for
imaginary part of z^2 equal to 0.0 or -0.0, if the Lisp implementation recognizes
signed zero.
Resolution of #1 and #3 is to work in upper half plane. Resolution of #2 is to
replace (EXPT Z 2.0) with (* Z Z).
Incidentally need to increase tolerance for erf(-0.75+%i) test as (* Z Z) yields
slightly different results for some Lisps.
diff --git a/src/gamma.lisp b/src/gamma.lisp
index b56ec3d..3cbc7df 100755
--- a/src/gamma.lisp
+++ b/src/gamma.lisp
@@ -1940,7 +1940,7 @@
;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
;;
;; When z is real sqrt(z^2)/z is signum(z). For complex z,
-;; sqrt(z^2)/z = 1 if abs(arg(z)) <= %pi/2 and -1 otherwise.
+;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
;;
;; This relationship has serious round-off issues when z is small
;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
@@ -1949,14 +1949,26 @@
;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
;; taken to return real results for real arguments and imaginary
;; results for imaginary arguments
+;;
+;; Pure imaginary z with Im(z) < 0 causes trouble for Lisp implementations
+;; which recognize signed zero, so just avoid Im(z) < 0 altogether.
+
(defun complex-erf (z)
+ (if (< (imagpart z) 0.0)
+ (conjugate (complex-erf-upper-half-plane (conjugate z)))
+ (complex-erf-upper-half-plane z)))
+
+(defun complex-erf-upper-half-plane (z)
(let ((result
(*
- (if (< (realpart z) 0.0)
+ (if (< (realpart z) 0.0) ;; only test needed in upper half plane
-1
1)
(- 1.0
- (* (/ (sqrt (float pi))) (gamma-incomplete 0.5 (expt z 2.0)))))))
+ ;; GAMMA-INCOMPLETE returns conjugate when z is pure imaginary
+ ;; with Im(z) < 0 and Lisp implementation recognizes signed zero.
+ ;; Good thing we are in the upper half plane.
+ (* (/ (sqrt (float pi))) (gamma-incomplete 0.5 (* z z)))))))
(cond
((= (imagpart z) 0.0)
;; Pure real argument, the result is real
diff --git a/tests/rtest_gamma.mac b/tests/rtest_gamma.mac
index 3a9d910..f70975a 100755
--- a/tests/rtest_gamma.mac
+++ b/tests/rtest_gamma.mac
@@ -2490,7 +2490,7 @@ closeto(
erf(-0.75+%i),
(-1.372897192365736489613456241111589390954675856186764729607156305b0
+ 0.539788632227100129936591912063260716699852732091113612337142798b0*%i),
- 5.0e-15);
+ 6.0e-15);
true;
closeto(