gmp_printf bug?

Torbjorn Granlund tg at gmplib.org
Wed Aug 3 23:06:59 CEST 2011


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

  > I hacked my code to avoid the exp calls.  I suppose this is what you
  > suggested:
  
  Something like that. Looks nice and simple when you implement it.
  
Excep now it is not very obvious how it works.  :-)

  > Unfortunately, it stillruns stupidly slow at about 40 seconds on the
  > main GMP machine, using dumbmp.  Some improvements to dumbmp are needed.
  > With a less silly dumbmp mpz_root starting approximaton, time is down to
  > 2.3s, which is OK.
  
  So it is mpz_root which is the bottleneck (in the case of linking with
  dumpmp)? What is the less-silly approximation? nthroot(x, n) \appr
  2^{bitsize(x) / n} or so?
  
Yep, I find the smallest integer t, for which 2^t > nthroot(x,n).

  After a quick look at dumpmp.c, I guess there are also significant
  savings by having special cade to handle the case of *square* roots
  (call to mpz_pow_ui unnecessary, and call to mpz_tdiv_q_ui can be
  replaced by a shift).

I pushed such improvements an hour ago...
  
  BTW, do you have any error analysis? To me, it seems that if you want to
  guarantee that the returned value can never be larger than the true
  value of log_b(2) = log(2) / log(b), the square roots used, as well as
  the multiplications, should be rounded *upwards* to the given precision.
  Naturally, this doesn't matter much when the function is used only for a
  small number of inputs, and you have some independent check that the
  results are correct for all these inputs.
  
Indeed, the current code is not 100% right in this sense.

Given the limited use of the code, I was thinking of computing with many
extra bits, simply checking that no final value is close to a tie case,
and abort if that happens.

The code should also generate floor(B*log(b)/log(2)).  Then we can get
rid of all floating point arithmetic in GMP, including very slow
floating point divides.

The other main purpose is to make printing of huge mpf values work
properly.

-- 
Torbjörn


More information about the gmp-devel mailing list