Anomaly in mpn_sqrtrem and mpn_rootrem
tg at gmplib.org
Wed Jul 15 10:46:19 UTC 2015
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.
Please encrypt, key id 0xC8601622
More information about the gmp-devel