gmp_printf bug?

[Zimmermann Paul]

       Torbjörn,

> The problems is in mpf/get_str.c.  There is the line,
> 
> e = (unsigned long) n_more_limbs_needed * (GMP_NUMB_BITS * mp_bases[base].chars_per_bit_exactly);
> 
> which computes an exponent for use by mpn_pow_1_highpart.  I am looking
> at the 20 lines of code staring around line 180.  The value in
> n_more_limbs_needed is large, with the bug report exponent, it becomes
> 21022920054702079.
> 
> With the approximation of log(2)/log(10) in chars_per_bit_exactly, e
> becomes 405025890106316032.  If it had been computed exactly, it would
> have become 405025890106316169, i.e., 137 greater.
> 
> Any suggestions?  Calling mpf recursively to do the floating point
> arithmetic more accurately will not be considered a realistic solution.
> :-)

one solution would be to replace chars_per_bit_exactly by a rational
approximation p/q, with p and q two 32-bit integers. Then you can get
an accuracy of about 64 bits on p/q.

For example log(2)/log(10) is well approximated by 579001193/1923400330,
with error less than 10^(-20).

Then you can estimate e with a call to umul_ppmm and a division of two limbs
by one limb.

If you want to avoid the division, you can use a 64-bit value for p, and
a power of two for q.

Paul