nisse at lysator.liu.se
Thu Sep 5 14:32:38 UTC 2019
tg at 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
> 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
I think remainder needs to be principal (could possibly allow equality
with the divisor), since we alternate between
a <-- a - q b
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
a -= b;
if (a is too small), exit
if (a >= b)
q = a / b;
a -= q * b;
if (a too small)
a += q; q--;
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.
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
More information about the gmp-devel