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

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

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/