bug in addrr

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

I've just noticed the following bug in addrr() [ has been with us since the
beginning ]:

(16:28) gp > \p28
   realprecision = 28 significant digits
(16:28) gp > -1. + precision(0.99999999, 3)
%1 = -0.00000001001171767711639404
(16:30) gp > \x
[&=08332864] REAL(lg=4,CLONE):05000004  (-,expo=-27):c07fffe5  ac000000 00000000
(16:31) gp > precision(%1,3)
%2 = -0.0000000100117176

It is bad enough that PARI's accuracy is necessarily a multiple of 32bits,
so that we can't do better than %2 (we have one correct bit, but do as if we
had 32).

But %1 is even worse: an extra spurious word is added.

The real fun starts when you repeat this kind of things [at \p28]:

  x = (1. + 10^-28) - 1;
  for (j=1, 100, x = (x-10^-28) + 10^-28)

at the end we still have x ~ 1e-28 (correct), but x boasts 201 words of
precision (of which about a single bit can be trusted).

I've fixed this in CVS [can't fix it in mp.s: arch 68k still has the bug...]
The fix is not perfect: we lose a word of precision now and then (but most
of the time, only the first 2 or 3 leading bits were correct). The only clean
way out I see is to add a component to the t_REAL type to specify the exact
relative precision (in bits, this time).

=============================================================================

Beside the fact that precision is given in words (using bits would solve
everything), the problem is due to the fact that PARI insists on floating
point numbers (t_REAL) being "normalized", i.e that the highest word in the
mantissa has its highest bit set to 1 [ if necessary, the limbs are shifted,
and exponent updated ]

I think the following representation would be nicer (and be easier to code):
drop the condition that mantissa be normalized and impose instead that the
exponent of any t_REAL is a multiple of BITS_IN_LONG [ better: store exponent
with respect to basis 2^BITS_IN_LONG ]  (I believe that is basically GMP's
choice)

I see three places where the current representation is better:

1) multiplication by any power of 2 is trivial ( esp. in place: setexpo() )

2) computing the exponent of a t_REAL is trivial [simpler than for a t_INT,
see expi(), where bfffo() has to be called].

3) division doesn't need a preliminary normalization (shift) of the divisor


On the other hand, currently:

1) t_REAL are constantly shifted back and forth (on average 0.5 shifts per
multiplication, 1.25 shifts per addition, 1.5 shifts per subtraction)

2) The code in addrr() becomes quite involved (read: inefficient) to prevent
uncontrolable precision increase like the above. A spurious new word of
precision is added each time we shift even a single bit (provided the shift
is not a multiple of the word length). We have to remember this, and possibly
compensate before returning (if final normalizing shift offsets the initial
one). For short inputs (say 2 or 3 words), shifts (and auxiliary code) might
even dominate the addition time.


In the PARI code, additions/subtraction are much more frequent than in place
multiplication by powers of 2. As for expo(): the "trivial" exponent would be
in base 2^BITS_IN_LONG, and the (more precise) exponent currently store would
need to be recovered using expi(). I think most crucial calls (certainly all
the ones in mp.c) could be updated to use an exponent in base 2^BITS_IN_LONG.

Before I seriously consider changing this representation, does anybody have
an idea where else (or what for) it could be useful ?

    Karim.

P.S: The change would be relatively easy to do if we dump mp.s (which would
be a good thing anyway since we cannot test it anymore. The generic kernel
should be good enough... Also it would simplify integrating a GMP kernel,
which would then provide even better 68k support:-)
-- 
Karim Belabas                    Tel: (+33) (0)1 69 15 57 48
Dép. de Mathematiques, Bat. 425  Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud             Email: Karim.Belabas@math.u-psud.fr
F-91405 Orsay (France)           http://www.math.u-psud.fr/~belabas/
--
PARI/GP Home Page: http://www.parigp-home.de/