Karim BELABAS on Fri, 4 Oct 2002 15:33:15 +0200 (MEST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: bug in pari-gp precision? (fwd) |
Hi, I'm forwarding a very interesting message I received about floating point accuracy in PARI. Two quick comments: 1) the problems reported are not due to my addrr() patches: I've checked version 1.39.15, they were already there. 2) no worries about polroots: computations and checks are done using exact arithmetic. 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/ ---------- Forwarded message ---------- From: Walter Neumann <neumann@cpw.math.columbia.edu> Date: Fri, 4 Oct 2002 02:21:32 -0400 (EDT) To: Karim BELABAS <Karim.Belabas@math.u-psud.fr> Dear Karim A new one: ? 1. + 10^-50 - 1. %2 = -2.52435489 E-29 What is happening here? I visited pari-dev's archive and see there has been no discussion of your post on the precision issue. You ask there if a correction would have application. It would be very important for us. Let me first say how we've been using pari: Oliver Goodman, Craig Hodgson, and I developed, over the last eight years, a 3-manifold package called Snap, based on Jeff Week's Snappea, using the pari libraries, to do exact provable computations for hyperbolic 3-manifolds, computing hyperbolic structures and arithmetic invariants of them. We promise (in a paper in Experimental Mathematics vol 9 (2000)) that Snap gives proofs of the hyperbolic structures it finds. The "proof" involves giving an exact ideal hyperbolic triangulation of the manifold in question. Each ideal simplex is given by a complex parameter which is an algebraic number. We thus describe the parameter by giving an an approximate numerical value (at least 50 digits pari precision by default) plus an exact description as an element of an appropriate number field. As I understand pari's real arithmetic, although we can be supremely confident of Snap's results, our claim of "proof" is wrong, and is hard to correct with current pari. The numerical value of the parameter is computed as a root of an integer polynomial and we write: 'we quote from the manual for the pari libraries: ``The algorithm used is a modification of A. Sch\"onhage's remarkable root-finding algorithm, due to and implemented by X. Gourdon. Barring bugs, it is guaranteed to converge and to give the roots to the desired accuracy.'' This guarantee seems to me suspect? Moreover, for geometric reasons we need to verify that the parameters have positive imaginary part, which we do from the numerical value. We must also confirm that the equations that say the tetrahedra fit together correctly are exactly satisfied. Numerical issues come up here because, these equations say that certain products of the parameters equal 1 (easily verified exactly) but also that the sum of their logs is 2\pi*i which we verify numerically (the precision need only be better than \pi). But, without analyzing carefully a complicated numerical computation that pari has done, one has no proof that pari's answer and claimed precision are even close to the truth, as in pari's ? (1. + 10.^-20 - 1. - 10.^-20 + 10.^-27 - 10.^-27 )*10.^29 + 1. %13 = -0.2621774483536188886 ? \x [&=080e1eec] REAL(lg=4,CLONE):05000004 (-,expo=-2):c07ffffe 863c1f5c dae42f94 We really need both the number theoretic features of pari and the real arithmetic (the same is true for the undergraduate student Max Lipyanskiy, who first raised this issue with you, in connection with an accurate Dirichlet domain program he has written). Since one of the points of our programs is that they are supposed to provide "proven" results, we need honest report of the valid precision of real arithmetic. Thus, I would welcome your suggestion of adding a bit-count of precision to the data-structure for a pari real, and having this bit-count treated honestly. Since it would be too conservative for "practical use" one would probably want to make it accessible to the program, but take a more pragmatic approach, as in the current implementation, to what one displays. Regards, Walter