Karim Belabas on Fri, 10 Jun 2005 14:53:03 +0200

 Re: precision question

```* Rampal Etienne [2005-06-10 14:09]:
> I have a problem with precision in PARI.
> Suppose we have default(realprecision,30).
> Then PARI gives
>
> > 1E+0 - 1E+0
> 0.E-38
>
> Is it possible to get
>
> > 1E+0 - 1E+0
> 0
>

No. See section 1.2.6.3 and 1.4 in the "User's guide to PARI/GP" for an
elaboration.

> In other words, if I want to calculate 1E+0 - 1E+0 + 1E-40, I want to get
> 1E-40 and not 0.E-38 which is the answer PARI gives me right now.
> This is important, because in a PARI program it could occur that the
> difference of exactly the same number is taken (e.g. log(10/2) - log(2+3) +
> 1E-40), which then gives the wrong answer.

It's not a wrong answer, it's a correct answer (= some number less than 1e-38),
which is not as precise as you might expect from an algebraic computation
involving infinite accuracy.

It's a common misconception about what floating point computations should
achieve. Floating point arithmetic is a precise model of computation for
finite number systems, within which one computes much faster than with
exact algebraic expressions, at the expense of weaker results (due to
rounding errors). Basic mathematical properties of real numbers do not
hold in floating point arithmetic: addition is not associative, a sequence
of additions is not commutative, etc.

A floating point number such as 1E+0 is not just a notational device for
an exact rational number writen in decimal (actually an integer here), it
represents _any_ real number not too far from 1. There is no reason why two
incarnations of 1E+0 should represent the same real number.

What you want [ log(10/2) - log(2+5) is an exact 0 ] can only be achieved by
a true computer algebra system (which PARI is not), not by any floating
point calculation.

Hope this helps,

Karim.
--
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dep. de Mathematiques, Bat. 425   Fax: (+33) (0)1 69 15 60 19
Universite Paris-Sud              http://www.math.u-psud.fr/~belabas/
F-91405 Orsay (France)            http://pari.math.u-bordeaux.fr/  [PARI/GP]

```