Karim.Belabas on Wed, 4 Apr 2001 11:16:41 +0200 (MET DST)

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: [GP/PARI-2.1.0] arithmetic weirdness

On Wed, 4 Apr 2001, Gerhard Niklasch wrote:
> 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.

This is more or less what is being done internally. All transcendental
functions are computed with one extra word of precision, then the result is
truncated. The first two things (t_REAL) LLL does are
1) compute the number of significant words used in the input
2) add one guard word to all entries [ x = gmul(x, realun(prec+1)) ]

It's a(n arguable) design choice, for the sake of maximal speed for low
precision inputs: the t_REAL type corresponds to low-level multiplication;
correct rounding and precision settings are the caller's responsibility.

Of course you may end up with a number of recursive calls adding up their own
extra word of precision until sluggish death occurs [ remember the zetak
soaring accuracy problem ? ], so this places another burden on the PARI
programmer [ stack management being the number 1 annoyance, until you get
really used to it. Ubiquitous typecasts as opposed to, eg. structs or OO
approach being my number 2 ].

It would be nicer if PARI also supported directed rounding or interval
arithmetic; very useful in a number of settings. Such a package was submitted
3 years ago by Manfred Radimersky but I don't like changing the kernel [e.g
Karatsuba multiplication for t_REAL has been there for 4 years, and is still
not enabled by default, nor is it documented!] and, to my great regret, I
unfortunately never had to time to check and adapt/include it. It still is
on my personal TODO list.


Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
PARI/GP Home Page: http://www.parigp-home.de/