t...@gmplib.org (Torbjörn Granlund) writes:
I suppose a "naive" implementation scales the full dividend and the full
divisor by a factor which fits in a limb. One could in theory save some
time by scaling just their O(log(Q)) most significant limb.
That paragraph was supposed to read:
I suppose a "naive" Svoboda implementation scales the full dividend and
the full divisor by a factor which fits in a limb. One could in theory
save some time by scaling just their O(log(Q)) most significant limbs.
How likely does Svoboda overesimate a quotient limb?
Well, one could of course adjust down it in O(1) time by involving
another set of high limbs.
Here is a somewhat less broken version of my algorithm:
udiv_qr_3by2 (qnext, n1, n0, np[0], np[-1], np[-2], d1, d0, dinv);
for (...)
{
// Apply q to most significant partial remainder limbs.
// Think: Is size = 2 enough or do we need size = 3?
// Optimisation: Results from last divappr should already have done most of
// the computation here; we presumably need just one more mul+subtract.
cy = submul_1 (np - 2, dp + dn - 2, qnext, 2);
ASSERT (np[0] >= cy); /* FIXME: Possibly handle spillage */
np[0] -= cy;
q = qnext;
// Now, with the most significant few limbs of the next partial remainder
// computed to some accuracy, compute q for the next iteration. This will
// overestimate q somewhat more often than with the old algorithms.
udiv_qr_3by2 (qnext, n1, n0, np[0], np[-1], np[-2], d1, d0, dinv);
// Apply q to all but the most significant partial remainder limbs.
cy = submul_1 (np - dn, dp, q, dn - 2);
if (UNLIKELY (cy > np[-2]))
{
np[-2] -= cy;
mpn_sub_1 (np - 1, np - 1, 2, 1);
// q was too large, qnext is really off. Recompute it!
// (We could perhaps recompute it more cleverly.)
ASSERT_CARRY(mpn_add_n (np - dn, np - dn, dp, dn));
udiv_qr_3by2 (qnext, n1, n0, np[0], np[-1], np[-2], d1, d0, dinv);
}
else
np[-2] -= cy;
np--;
}
--
Torbjörn
Please encrypt, key id 0xC8601622
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
