Re: [GP/PARI-2.1.0] arithmetic weirdness

[Gerhard Niklasch <nikl@mathematik.tu-muenchen.de>]

P.S. to
> Message-Id: <200104032037.WAA02692@sunscheurle3.mathematik.tu-muenchen.de>
> From: Gerhard Niklasch <nikl@mathematik.tu-muenchen.de>
> Date: Tue, 3 Apr 2001 22:37:57 +0200 (MEST)

> (22:16) gp > 1/100.
> %49 = 0.009999999999999999999999999999
> (22:16) gp > \x
> [&=004a9bcc] REAL(lg=5,CLONE):05000005  (+,expo=-7):407ffff9  a3d70a3d  70a3d70a  3d70a3d7  
> 
> (note that this is *correctly* rounded -- the next four bits would
> all be zero!)
...
> (22:13) gp > ((2^144\/100)/2^144+0.)
> %47 = 0.009999999999999999999999999999
> (22:15) gp > \x
> [&=004a9b7c] REAL(lg=5,CLONE):05000005  (+,expo=-7):407ffff9  a3d70a3d  70a3d70a  3d70a3d7  
> 
> (correctly rounded, but _looks_ no better than the original)
...
> and, moreover,
> 
> (22:22) gp > %47+1/2^102
> %55 = 0.01000000000000000000000000000
> (22:22) gp > \x
> [&=004a9cbc] REAL(lg=5,CLONE):05000005  (+,expo=-7):407ffff9  a3d70a3d  70a3d70a  3d70a3d8  
> (22:22) gp > %55*100
> %56 = 1.000000000000000000000000000
> (22:22) gp > \x
> [&=004a9ce4] REAL(lg=5,CLONE):05000005  (+,expo=0):40800000  80000000  00000000  00000000  
> (22:23) gp > %47*100 
> %57 = 0.9999999999999999999999999999
> (22:24) gp > \x
> [&=004a9d0c] REAL(lg=5,CLONE):05000005  (+,expo=-1):407fffff  ffffffff  ffffffff  ffffffff  
> 
> i.e. with *correct* rounding at the default realprecision, 100*(1/100.)
> does not numerically equal 1+0. ...

I shouldn't be posting late at night, and this is *not* correctly
rounded.  Ugh.

The exact result of the multiplication %47*100 would have five extra
1 bits to the right, which we're pushing out during normalization.

We shouldn't drop them, we should remember enough of them to round
upward.  The standard guard bits + sticky bit technique is being
called for here.

For elementary operations and for as many algebraic and transcen-
dental functions as feasible, we should return the correctly rounded 
t_REAL representation, at the appropriate precision, of the exact
result of applying the operation  (in the reals, not in the t_REALs)
to the input value(s).  Perhaps stuff like LLL would be a trifle
better behaved if we didn't drop bits in simple multiplications...
It is not necessary to compute the full 192-bit product of two 96-bit
significands to achieve this;  a few extra bits suffice to do the
tracking.

Cheers, Gerhard  (who gave away his old numerical analysis textbooks
to students years ago - but, as I discovered whilst exchanging a couple
mails with Will Galway off the list last night, typing `IEEE754' at
www.google.com pulls up a lot of useful material)