gmp_printf bug?

[Torbjorn Granlund]

nisse at lysator.liu.se (Niels Möller) writes:

  Torbjorn Granlund <tg at gmplib.org> writes:
  
  > I am not sure how to compute the more-than 53 bits accurate logarithm
  > tables needed for solution 1, in the very constrained environment of a
  > GMP compile.
  
  Is it possible to compute the 53-bit logarithm and improve it
  iteratively? I don't know the standard way of computing logarithms
  numerically; since its a transcendental function, I guess there's no
  convenient newton iteration (but if exp is easier to compute, there
  ought to be an iteration using exp evaluation and newton).
  
We know well how to compute exp with integer base and exponent, but it
is not as easy when the exponent is a fraction, as it will be here.

Since we only need a few dozen bits, we don't need quadratic
convergence; adding a bit at a time will work.  We just need
exponentiation.

Fortunately, we need exponentiation in a quite restricted sense.  We
just need to compute b^e for integers b, 3 <= b <= 255 where the
expected result is close to 2, i.e., e is an approximate b-logarithm of
2.

We may scale e to become an integer by multiplying it by a power of 2,
and then perform exponentiation.  The result will be a forbiddingly huge
number, but we're really only interested in the high word(s).  We can
therefore light-heartedly throw away low words during exponentiation,
and just keep the most significant word plus some dozen bits.  (I write
words, not limbs, since these need to be size_t large quanities.)


I have convinced myself that we need to do away with the f-p table in
mp_bases for good.  We rely on the host system's log function's
exactness, which is fragile.  The numbers will be some ULP too small or
some ULP too large, we know not.  I mpf/get_str.c, we need to carefully
round up or down when using these precomputed values, the direction
depends on the exponent of the number to be converted to a string.

With integer (fractions) representing the logarithms, we get more
adequate accuracy, and we can make them right to the last bit.  We can
easily get rounding down (table the chopped log value, use umul_ppmm and
extract high word).  Rounding up will be slightly tricker, using the
chopped log value + 1 is not quite enough, if we chop the result product
implicitly by taking the upper umul_ppmm word.

-- 
Torbjörn