float evaluation of erfi



>>>>> "Robert" == Robert Dodier <robert.dodier at gmail.com> writes:

    Robert> 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. 

    Robert> Without bringing signed zero into the picture, can you show that
    Robert> sqrt((-i)^2)/(-i) is something other than -1 ?

Without signed zeroes, of course, -1 is the only answer.

Some of the numerical algorithms will probably need to be examined
more closely.  One of the first things the numerical algorithms do is
convert x+y*%i into the equivalent of (complex x y).  But maxima has
pure imaginary numbers y*%i, but CL does not, especially if your lisp
has signed zeroes.

    Robert> Given erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2, z^2)/sqrt(pi)),
    Robert> if sqrt(z^2)/z is 1 for pure imaginary z, then erf(-z) = erf(z) (which
    Robert> is how I bumped into this to begin with). I don't think we want that.

Agreed.  We should take more advantage of known properties in the
numerical routines like erf(conjugate(z)) = conjugate(erf(z)) and
erf(-z) = -erf(z).  Then the numerical routine can just concentrate on
getting the correct value of erf in the first quadrant.

Ray


    >> CL doesn't forbid signed zeroes.

    Robert> Right, but it doesn't require them, so any code depending on them is
    Robert> nonportable.

    Robert> For the record, I've committed and pushed 8f10856b70 to address the
    Robert> problems with erf & friends.

I have not yet checked, but does a similar change need to be made for
bfloats?  Unfortunate that there's a slight loss of accuracy with this
change.  Oh well.

Ray