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/