Karim BELABAS on Thu, 26 Sep 2002 01:06:10 +0200 (MEST)

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

bug in addrr

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

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 ?


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/