Anomaly in mpn_sqrtrem and mpn_rootrem

[Torbjörn Granlund]

bodrato at mail.dm.unipi.it 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.


-- 
Torbjörn
Please encrypt, key id 0xC8601622