Gerhard Niklasch on Wed, 4 Apr 2001 07:31:52 +0200 (MEST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: [GP/PARI-2.1.0] arithmetic weirdness |
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)