On 3/18/09, Raymond Toy wrote:
-------------------------------------------------
Robert Dodier wrote:
> On 3/17/09, Raymond Toy <toy.raymond at gmail.com> wrote:
>
>
>> I can reproduce this problem with current CVS using gcl 2.6.8 on
>> Solaris. It's a bit slow on g(8,40), but it gets some kind of
>> unrecoverable error (fault count too high) with g(9,40) and proceeds to
>> dump core.
>>
>> This doesn't happen with cmucl. I don't know why gcl would fail in this
>> way.
>>
>
> Looks like exp(-1/(1 - xx)) underflows in an ungraceful way ....
> g(7, 40) yields a term exp(- 9.999999999712443b6) =>
> 1.51737304873662b-4342945, which (I guess) causes trouble.
> With Clisp it triggers the error "ASH: too large shift amount 14426900".
> The error is in FPPLUS. Is that the bigfloat code is trying to
> normalize terms to a common basis or something? Just a guess.
>
>
Interesting. My version of clisp (2.47) goes ahead and computes the
result. And a different version of gcl (2.6.8 pre, whatever that
means), works just fine too.
Ray
-----------------------------------------------------------------
I have checked that this same bfloat underflow hang occurs with the
same code in versions 5.9.3, 5.12.0, 5.13.0, 5.14.0, 5.15.0, 5.16.1, as well
as
5.17.3 (all GCL 2.6.8 except for v. 5.9.3 which uses GCL 2.6.7).
I suppose a workaround is to define my own exponential function with
graceful
underflow behavior, but isn't this something which should be part of the
bfloat
code??
Ted Woollett
-------------------------------------------------------------
tests done with batch file bfbug.mac:
------------
display2d:false$
h(xx) := exp(-(xx+1)/4 )/4 + exp( -1/(1-xx))/(1-xx)^2$
" test h(xx) with default floats "$
gg(j ) := block([x,y],
y : float(10^(-j)),
x : float(1-y),
print(j,y,x),
float(h(x)) )$
fpprec;
gg(7);
gg(8);
" test h(xx) with bfloat and fpprec:40 locally "$
g(j,fp) := block([fpprec,x,y],
fpprec:fp,
y : bfloat(10^(-j)),
x : bfloat(1-y),
print(j,y,x),
bfloat(h(x)) )$
g(7,40);
g(8,40);
------------------------
Maxima 5.17.1 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (aka GCL)
(%i1) batch(bfbug);
batching #pc:/work3/bfbug.mac
(%i2) display2d : false
(%i3) h(xx):=exp((-1)/(1-xx))/(1-xx)^2+exp((-(1+xx))/4)/4
(%i4) " test h(xx) with default floats "
(%i5) gg(j):=block([x,y],y:float(1/10^j),x:float(1-y),print(j,y,x),
float(h(x)))
(%i6) fpprec
(%o6) 16
(%i7) gg(7)
7 9.9999999999999995E-8 0.9999999
(%o7) 0.15163266871898
(%i8) gg(8)
8 1.0E-8 0.99999999
(%o8) 0.15163266530724
(%i9) " test h(xx) with bfloat and fpprec:40 locally "
(%i10) g(j,fp):=block([fpprec,x,y],fpprec:fp,y:bfloat(1/10^j),x:bfloat(1-y),
print(j,y,x),bfloat(h(x)))
(%i11) g(7,40)
7 1.0b-7 9.999999b-1
(%o11) 1.516326687189750264901169661977624778275b-1
(%i12) g(8,40)
8 1.0b-8 9.9999999b-1
(Maxima hangs and must restart )
-------------------------------