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 ]


[ 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/