Re: [GP/PARI-2.1.0] arithmetic weirdness

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

In response to:
> Message-ID: <20010403220706.U7065@yellowpig>
> Date: Tue, 3 Apr 2001 22:07:06 +0200
> To: Philippe Elbaz-Vincent <pev@pev.math.univ-montp2.fr>
> Cc: PARI developers <pari-dev@list.cr.yp.to>
> References: <200104020600.IAA01232@sunscheurle3.mathematik.tu-muenchen.de> <Pine.LNX.4.33.0104020839360.6239-100000@pev.math.univ-montp2.fr>
> From: Bill Allombert <allomber@math.u-bordeaux.fr>
> 
> I feel it is indeed a bug in PARI :
> 
> consider the following computation:
> 
> ? (2^96\10)/2^96+0.
> %10 = 0.09999999999999999999999999999
> ? (2^96\/10)/2^96+0.
> %11 = 0.1000000000000000000000000000
> 
> The second result seem more accurate, since we
> use the smallest residue.
...
> (the 0x9 has been rounded up to 0xa, which is better since next hexa is a 9)

Good observation.  I missed that.  Ok, this should be fixed then.

However, while this will give the expected (for some values of
expect :) result here, surprises elsewhere are still unavoidable,
for the fundamental reason that "most" rationals aren't exactly
representable.  Users will still need to be warned of this  (errm,
they *are* being warned of this, or were last time *I* read the
documentation  (long ago)).  E.g.

(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:15) gp > ((2^96\/100)/2^96+0.)
%48 = 0.009999999999999999999999999995
(22:16) gp > \x
[&=004a9ba4] REAL(lg=5,CLONE):05000005  (+,expo=-7):407ffff9  a3d70a3d  70a3d70a  3d70a3c0  

(worse -- the proposed fix isn't quite right, it doesn't round at the
correct bit position.  lgefint is beside the point;  you need the number
of bits occupied by 10^n, and the number of words implied by the current
realprecision)

(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)

(NB the binary->decimal conversion might also warrant another look)

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. , and the representable x which comes
closest to solving 100*x==1+0. numerically is not numerically equal to
the correctly rounded binary result of the division 1/100. at our default
precision.

Such is life in binary.

(Or in any other finite base, to finite precision.)

Cheers, Gerhard