>>>>> "Robert" == Robert Dodier <robert.dodier at gmail.com> writes:
Robert> There is some strangeness in COMPLEX-ERF
Robert> (src/gamma.lisp). I hope I'm not mistaken about it. Oh
Robert> well, here goes.
Robert> (1) This comment
Robert> ;; When z is real sqrt(z^2)/z is signum(z). For complex
Robert> ;; z, sqrt(z^2)/z = 1 if abs(arg(z)) <= %pi/2 and -1
Robert> ;; otherwise.
Robert> is incorrect, right? When z is pure imaginary and Im(z) <
Robert> 0, sqrt(z^2)/z is equal to -1, right? (sqrt((a i)^2)/(a i)
Robert> = sqrt(a^2)/a = signum(a), right?) So i think the sign
Robert> computation
This is my reasoning. For a > 0, let z = a*i = a*exp(%i*%pi/2). Then
z^2 = a^2*exp(%i*%pi) so sqrt(z^2) = a*exp(%i*%pi/2) = z. Thus
sqrt(z^2)/z = 1. For a < 0, we have z = |a|*exp(-%i*%pi/2), so
sqrt(z^2) = sqrt(|a|^2*exp(+%i*%pi)) = |a|*exp(+%i*%pi/2). Thus
sqrt(z^2) = |a|*exp(+%i*%pi/2)/(a*exp(-%i*%pi/2)) = |a|/a *exp(%i*%pi)
= 1.
However, you can quibble over whether exp(-%i*%pi) is exp(%i*%pi)
since that's on the negative axis.
I cannot find it right now, but there is some relationship on
functions.wolfram.com that mentioned f(z) = sqrt(z^2)/z*g(z) and also
said f(z) = g(z) if realpart(z) >= 0 and f(z) = -g(z) if realpart(z) <
0.
Robert> (2) Another problem.
Robert> (gamma-incomplete 0.5 (expt z 2.0))
Robert> Some Lisp implmentations (CMUCL, SBCL, ECL, GCL, not
Robert> Clisp) compute (EXPT Z 2.0) as having a nonzero imaginary
Robert> part when Z is pure imaginary. That causes trouble for
For cmucl, I think this is a bug. There are compiler transforms so
that (expt z 2) is computed using (* z z). But from the repl (expt z
2) uses logs. That's quite annoying and should be fixed.
Robert> GAMMA-INCOMPLETE because, depending on the sign of the
Robert> spurious bit, the return value might be the conjugate of
Robert> the correct value (which you get for imaginary part = 0.)
Robert> It just so happens that for all the Lisp implementations
Robert> which have the spurious bit, it has a sign such that the
Robert> result comes back as the conjugate. But then we multiply
Robert> the imaginary part by the wrong sign -- et voila, two
Robert> wrongs do make a right!
Heh. Getting branch cuts correct is very difficult. Not having a
signed zero makes it even harder.
Anyone interested in this should read Kahan's article Much Ado About
Nothing's Sign. (Some people disagree with Kahan, though.)
Robert> Not sure what to do at this point. Any code which depends
Robert> on negative zero is nonportable, as negative zero doesn't
Robert> exist in Common Lisp's model of floating point numbers. So
Robert> if there is some flag that we can use to suppress it, we
Robert> could enable that. Or in this limited context we could
CMUCL has no such flag. And since IEEE754 requires signed-zeroes it,
the only way to work around it is to add 0 to the result. (Because
-0.0 + 0.0 is required to be +0.0.)
Robert> Incidentally you can see what's going on by tracing
Robert> SIMP-ERFI, COMPLEX-ERF and GAMMA-INCOMPLETE, then
Robert> evaluating erfi(-0.15) and erfi(0.15).
I have a much simpler solution to erfi(-0.15). erfi is an odd
function, so compute erfi(-z) as -erfi(z).
That doesn't solve the general problem with branch cuts for
gamma-incomplete.
Ray