Paul.Zimmermann at loria.fr
Wed Jul 20 16:08:42 CEST 2011
> 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
> 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.
More information about the gmp-bugs