>>>>> "Robert" == Robert Dodier <robert.dodier at gmail.com> writes:
Robert> On 6/7/07, Barton Willis <willisb at unk.edu> wrote:
>> (%i65) ff:exp(-x^2/50^2)*sin(x)^2*cos(2*x)$
>> (%i66) gg : integrate(ff,x,-4*%pi,4*%pi)$
>> (%i67) float(gg);
>> (%o67) 0.0
>> (%i68) float(expand(gg));
>> (%o68) -6.153361529128077
>>
>> Look at both gg and expand(gg). You should see why (%o67) = 0.0, but
>> (%o68) isn't 0.0.
Robert> %o67 = 0.0 looks to me like a bug in GCL's floating point math.
Robert> There is one term exp(10000) and another one exp(-10000).
Robert> The first yields a floating point inf and the second, 0.0.
Robert> The product should be NaN (not-a-number).
Robert> SBCL could probably be coaxed into yielding NaN for %o67.
Robert> (Clisp cannot, to the best of my knowledge.)
:lisp (ext:set-floating-point-modes :traps nil)
for CMUCL will turn off the traps. This won't give NaN for %o67
because maxima's simplifier converts any product containing 0 to 0,
independent of what the other terms might be.
Robert> About this stuff,
>> (%i164) integrate(ff,x,-4*%pi,4*%pi);
>> (%o164) -(%e^(-36)*sqrt(%pi)*(3*erf((4*%pi+18*%i)/3)-6*%e^27*erf((4*%pi
>> +9*%i)/3)-6*%e^27*erf((4*%pi-9*%i)/3)+3*erf((4*%pi-18*%i)/3)+6*%
>> e^36*erf((4*%pi)/3)))/8
>>
>> (%i165) float(%);
>> (%o165) -5.1390589659105944*10^-17*(3.0*erf(0.33333333333333*(18.0*%i
>> +12.56637061435917))-3.192289443610792*10^+12*
>> erf(0.33333333333333*(9.0*%i+12.56637061435917))-3.192289443610792*10^
>> +12*erf(0.33333333333333*(12.56637061435917-9.0*%i))+3.0*
>> erf(0.33333333333333*(12.56637061435917-18.0*%i))+2.5867389201337716*10^
>> +16)
Robert> I wonder if Maxima could output something nicer here.
Robert> (1) Could erf applied to a complex float (i.e. float + %i * float)
Robert> yield a float or complex float?
I looked sometime ago for a routine to evaluate erf at complex
points. The only thing I found was something contributed to octave.
But I wasn't comfortable with it at the time.
Robert> (2) I tried abs(%o165) and realpart(%o165) and imagpart(%o165)
Robert> to try to get real floats from the %o165 but abs and friends didn't
Robert> seem to know how to do something useful with it. I wonder what
Robert> we could do to change that. Maybe (2) is moot if (1) is implemented.
Is it possible to express the real and imaginary parts of erf in terms
of other (simple) functions that maxima knows about in a useful way?
I notice functions.wolfram.com gives a few expressions, but they don't
seem particularly helpful.
Ray