| Karim BELABAS on Wed, 15 Nov 2000 12:53:53 +0100 (MET) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| Re: [GP/PARI] 2.0.21 bug in exp()? |
[Michael Somos:]
> the 2.0.21 and tried the following simple function :
>
> gp> f(x,d=10,v=0)=
> {
> local(rc);
> /*DEB*/ if(v,print("f("x","d")"));
> rc=
> if(x<=0, error("only positive allowed in f()"),
> if(d<0, error("recur depth positive"),
> if(x==1|d==0, 0,
> if(x<1, -f(1/x,d-1,v),
> /*x>1*/ exp(f(x-1,d-1,v))
> ))));
> /*DEB*/ if(v,print("f("x","d")="rc));
> rc;
> } /* end f() */
> gp> f(16/5,8,1)
[...]
> f(6/5,6)=4.289076188000423375581902986 E-1656521
[...]
> At this point the system seems to hang for a long long time. I can only
> guess that the bug may be in exp() somehow. Who knows? Shalom, Michael
At this point, PARI attempts to compute
exp(x = 4.289076188000423375581902986 E-1656521)
which is internally done as 1 + y, y = exp(x) - 1). Of course, y ~ x is
rather small. But the computation 1 + y can be done exactly since `1' is an
`exact 1' as opposed to `1.'
By PARI's philosophy that floating point operations always give a result
which is as precise as the input permits, 1 + y is a floating point number
whose mantissa has roughly 1656521 digits. Then things start to slow down...
It's quite difficult to correct this in a consistent and natural way. I could
use an arbitrary cutoff point (so as to decide when a precision increase will
be ridiculous), but I don't like this solution. For the time being, I'd
rather let the user correct his input, than change the PARI "philosophy".
For instance, what are we to do on input exp(1e-400) ??? [check the output
using \x, probably not what you want, but who knows ??? ]
A related problem:
(12:40) gp > exp(-1E10000000)
*** exponent overflow
[ also the error message is a bit surprising, but this is computed as
1/exp(-x), and the routine where the error occurs has no way to know about
this ]
Karim.
[ Side note: the PARI way of updating (increasing or decreasing) the
precision of objects, which are themselves computed to the largest precision
attainable given the input, is probably a bad idea. It can cause the
precision of transcendal operations to skyrocket, or immediately lose all the
digits in a floating point LLL run and abort, where an approximate answer
would have been just as good. It's tough to change this. ]
--
Karim Belabas email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France) Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://www.parigp-home.de/