Anomaly in mpn_sqrtrem and mpn_rootrem

Torbjörn Granlund tg at
Wed Jul 15 10:46:19 UTC 2015

bodrato at writes:

  > But I spotted "% k" and "/ k" there, and those are very expensive,
  > unless you table inverses of k, of course...
  There is already a single (there are two currently, but we can avoid one)
  division by k in the code, it is used to compute the first single-bit
  estimate. The trick I suggest consists in adding some fractional digits of
  log2(n), by table lookup, then use log2(n)/k to estimate a few of the
  highest digits of n^(1/k), not just one.
  I still have to analyse the maximal error, but I believe we can use the
  trick to skip some iterations.
I agree that your code is an improvement.

For all of GMP, I think we should consider alternatives to explicit
division (and modulus).  If the numerator is limited, we can always
compute an accurate quotient with a multiply, provided we have a
reasonable divisor inverse.

My division paranoia might not be motivated for all chips, though.
32-bit chips tend to have fast enough division hardware,, and recent
64-bit x86 CPUs have fast 64/64 division (but slow 128/64 division).
We cannot then beat the hardware with Newton.

Why bother?  Our sqrt(rem) and root(rem) will sometimes be used for just
a few limbs, and then division as part of the "bookkeeping" will take a
very large fraction of the time.

Please encrypt, key id 0xC8601622

More information about the gmp-devel mailing list