Gerhard Niklasch on Tue, 3 Apr 2001 22:37:57 +0200 (MEST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: [GP/PARI-2.1.0] arithmetic weirdness |
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