hgcd1/2

[Niels Möller]

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
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.