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