t...@gmplib.org (Torbjörn Granlund) writes:
> Can div2 be done sloppily, e.g., by extracting the top 64 bits of the
> dividend and the corresponding (up to) 64 bits of the divisor, and
> invoke DIV1 on these bits?
I think so. Result will be at most one off, except in the unlikely case
that q > d1, so one adjustment step should suffice on the likely code
path. I think I have code like that, but with a/b rather than div1
(a/b). I may even have posted it to the list.
And due to the stop condition of the double-limb loop in hgcd2, I don't
think any normalization shift is needed, just q = div1(ah, bh) should be
fine.
> OK, one needs to return a non-negative remainder, so a little bit of
> adjustment might be needed. But does the remainder need to be
> principal?
I think remainder needs to be principal (could possibly allow equality
with the divisor), since we alternate between
a <-- a - q b
and
b <-- b - q a
and code assumes that q >= 1 in each step and everything is
non-negative. Allowing an approximative (and too small) q is possible,
but then we need to add a conditional swap.
> Related to that: In hgcd2 you seem to handle a >= b just after a is set
> to a residue mod b. If that residue is principal, then clearly a >= b
> cannot happen. Perhaps I am reading the code wrong.
Maybe your confused by the code special casing q = 1? The iteration is
like
a -= b;
if (a is too small), exit
if (a >= b)
{
q = a / b;
a -= q * b;
if (a too small)
{
a += q; q--;
exit
}
}
at this point, we should always have a < b.
With a good enough div1, we should perhaps call it immediately, and get
rid of that branch too.
Regards,
/Niels
--
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
