There is some strangeness in COMPLEX-ERF (src/gamma.lisp). I hope I'm
not mistaken about it. Oh well, here goes.
(1) This comment
;; 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.
is incorrect, right? When z is pure imaginary and Im(z) < 0, sqrt(z^2)/z
is equal to -1, right? (sqrt((a i)^2)/(a i) = sqrt(a^2)/a = signum(a),
right?) So i think the sign computation
(if (< (realpart z) 0.0)
-1
1)
should be changed to
(if (or (< (realpart z) 0.0) (and (= (realpart z) 0.0) (< (imagpart z) 0.0)))
-1
1)
or some equivalent.
(2) Another problem.
(gamma-incomplete 0.5 (expt z 2.0))
Some Lisp implmentations (CMUCL, SBCL, ECL, GCL, not Clisp) compute
(EXPT Z 2.0) as having a nonzero imaginary part when Z is pure
imaginary. That causes trouble for GAMMA-INCOMPLETE because, depending
on the sign of the spurious bit, the return value might be the
conjugate of the correct value (which you get for imaginary part = 0.)
It just so happens that for all the Lisp implementations which have the
spurious bit, it has a sign such that the result comes back as the
conjugate. But then we multiply the imaginary part by the wrong sign --
et voila, two wrongs do make a right!
(3) One might think we could fix the sign problem by computing (* Z Z)
instead -- not so. It seems that the result has imaginary part equal to
negative 0.0 (the result of multiplying a negative number by 0.0) and
that is enough to change GAMMA-INCOMPLETE's return value.
Not sure what to do at this point. Any code which depends on negative
zero is nonportable, as negative zero doesn't exist in Common Lisp's
model of floating point numbers. So if there is some flag that we can
use to suppress it, we could enable that. Or in this limited context we
could invent a SAFE-SQR function to ensure the imaginary part is
unsigned zero. (EXPT Z 2) (i.e., integer exponent) doesn't help, btw.
Incidentally you can see what's going on by tracing SIMP-ERFI,
COMPLEX-ERF and GAMMA-INCOMPLETE, then evaluating erfi(-0.15) and
erfi(0.15).
Hope I haven't misunderstood everything. Apologies in advance if I have.
best
Robert Dodier