Re: [GP/PARI-2.1.0] arithmetic weirdness

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

In response to:
> Message-ID: <Pine.SOL.3.95.1010404105337.8039A-100000@geo>
> Date: Wed, 4 Apr 2001 11:16:41 +0200 (MET DST)
> From: "Karim.Belabas" <Karim.Belabas@math.u-psud.fr>
> To: pari-dev list <pari-dev@list.cr.yp.to>
> In-Reply-To: <200104040531.HAA02765@sunscheurle3.mathematik.tu-muenchen.de>
> 
> On Wed, 4 Apr 2001, Gerhard Niklasch wrote:
> [...]
> > 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).
> 
> 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.

Truncated?  (as opposed to rounded at the cutoff bit position?)

(there's many lines of code I've never looked at... :)

> 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.

What worries me is that t_REAL*t_REAL currently seems to be systemat-
ically biased  (towards zero),  which makes the caller's job a rather
difficult one, and the programmer's job  (predicting what accuracy is
going to be needed)  likewise.  The gp user or script programmer is
basically stuck with whatever gp calls to handle a * between t_REALs,
unless she resorts to rather devious tricks.

> 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 ? ],

Sure do. :)

Sticking on an extra word during each and every operation is not the
answer here -- we don't want to mess around all the time with numbers
which carry far more words around than warranted by the input precision,
and t_REAL wouldn't allow the recipient of a return value to see how
many of the bits are based on the input and how many are just padding.
The classic hardware FPU algorithm uses precisely three extra *bits*
(the last of which is sticky)  to achieve correct rounding for all
elementary operations, and the results returned are again rounded to
no longer than the input precision.

> It would be nicer if PARI also supported directed rounding or interval
> arithmetic; very useful in a number of settings.

Yes  (as would be many other things).  However, my current impression
is that we have neutral rounding for addition/subtraction and directed
rounding for multiplication, which is, well, a somewhat exotic com-
bination, and probably not helpful to iterative numerical algorithms
whether part of PARI or written as user scripts.

Cheers, Gerhard