Karim BELABAS on Sat, 5 Oct 2002 00:12:34 +0200 (MEST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: bug in pari-gp precision? (fwd) |
On Fri, 4 Oct 2002, Ilya Zakharevich wrote: > On Fri, Oct 04, 2002 at 03:33:15PM +0200, Karim BELABAS wrote: >> ? 1. + 10^-50 - 1. >> %2 = -2.52435489 E-29 >> >> What is happening here? > > Personally, I see no problem here (although I cannot guess what the > particular implementation is doing). You are running this with the > precision 28digits, right? PARI translates its "programs" from > strings to data "on the fly". It sees "1.", and creates a value with > 28 digits of precision. This means the absolute error is of order of > magnitude of 1e-28 (or 1-e29, depending on how you count ;-). If x > 0, it sounds sensible that x + y >= y for any y. Since PARI's floating point operations are not specified anywhere, it can't directly be called a bug. As you say, the relative error is correct. On the other hand, if you follow the computation, this is the result of "incorrect" rounding in one step, in the sense that there is a closest match than the one chosen, given the available data. Now, this goes against the philosophy of getting as precise answers as inputs warrant. So it is a bug in this respect. I am not 100% sure yet whether addrr(), mulrr(), divrr() should go to extra length to ensure proper rounding [ they probably should, but the cost would not be negligible at low accuracies ]. On the other hand, affir() and affrr(), at least, ought not to discard extra significant digits. After fixing these last two, the above phenomenon disappears. It is easy to fix divrr() also, without slowing it down. It suppresses most of the oddities when inputing floating point numbers. In particular the following [see the "arithmetic weirdness" thread from last year, starting pari-dev-1112]: (23:22) gp > \p38 (23:22) gp > 1e-20; (23:22) gp > \x [...] 88f4bb1c a6bcf584 (23:20) gp > \p28 (23:21) gp > 1e-20; (23:21) gp > \x [...] 88f4bb1c ^---- should round to d. I didn't notice any slowdown with these two changes. In any case, the kernel should be properly documented [ currently, the manual says things like: "affir(x,y): assigns the integer x into the real y", which is not helpful. ] I have added 3 TODO entries related to this: (Kernel section). 5 ensure consistent rounding in floating point operations. At least document it. Currently most routines truncate instead of rounding, most of the time. 4 support standard rounding modes in floating point operations 3 interval arithmetic It will probably be easier to port a GMP kernel instead. Cheers, Karim. -- Karim Belabas Tel: (+33) (0)1 69 15 57 48 Dép. de Mathematiques, Bat. 425 Fax: (+33) (0)1 69 15 60 19 Université Paris-Sud Email: Karim.Belabas@math.u-psud.fr F-91405 Orsay (France) http://www.math.u-psud.fr/~belabas/ -- PARI/GP Home Page: http://www.parigp-home.de/