Re: bug in pari-gp precision? (fwd)

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

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