Cofactor canonicalisation of mpn_gcdext

[Niels Möller]

Bill Allombert <Bill.Allombert at math.u-bordeaux1.fr> writes:

> 3) In my view there is a large use-case for mpn_gcdext where 0 <= s < |v|/g
> is suitable: when the user performs modular inversion, because it is assumed
> that g=1, and the value of t is not needed

I think that is why I have strived, when writing the new gcdext code,
to keep s >= 0. For some kind of symmetry, I have also tried to get t
<= 0, and this is why, in the case that v divides u, the solution s =
v, t = - (u - 1) is used rather than the "trivial" solution s = 0, t = 1.
That said, I don't object to the proposed 2|s| < |v|/g
canonicalization rule.

One other observation, before I get to the main part of this posting:
It is clear that there exists an s in the range 2 |s| < |v/g|, but it
wasn't clear to me that this is what's naturally produced by Euclid's
algorithm; I'm more familier with the bounds |s| < |v|/g, or the
stricter bound |s| + |some related matrix element| < |v|/g. If this
indeed holds, that means that when computing gcdext of two unsigned
single-word numbers, the cofactors fit nicely in signed single-word
numbers.

To get the old behaviour back, there are a couple of not too difficult
things that have to be done:

1. Change interface and implementation of gcdext_1, currently defined
   (for the odd case) as:

     /* Returns g, u and v such that g = u A - v B. There are three
        different cases for the result:
     
          g = u A - v B, 0 < u < b, 0 < v < a
          g = A          u = 1, v = 0
          g = B          u = B, v = A - 1
     
        We always return with 0 < u <= b, 0 <= v < a.
     */
     #if GCDEXT_1_USE_BINARY
     
     static mp_limb_t
     gcdext_1_odd (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b);

   in mpn/gcdext_1.c. This should be changed to returned signed
   cofactors using the desired canonicalization. Current code uses the
   binary algorithm here; I can think of some approaches to generate
   cofactors in the desired range using binary right-to-left
   processing, but the best way to do it is not obvious to me. So for
   the time being, I think it's best to revert to Euclid for this
   case. There's another version of gcdext_1, as the last function in
   that file (starting at line 271 and currently #if:ed out) which I
   think does the right thing, except for interface and the sign
   handling at return.

2. There's one call to gcdext_1, in mpn/gcdext_lehmer.c. This
   processing is the common case whenever the gcd is small, and it
   uses the cofactors returned by gcdext_1 to construct a linear
   combination of larger cofactors associated with the almost complete
   reduction (if that description makes any sense). Change this code
   adapt to the new gcdext_1 interface with signed cofactors.

   Since the absolute values of the small cofactors now fit in
   GMP_NUMB_BITS - 1 bits, one should use the new function
   mpn_addaddmul_1msb0 when available.
   
3. Hunt down the cases where gcdext returns early without a gcdext_1
   at the end of the reductions, and ensure that the return values
   conform to the canonicalization rules. I think the main place is in
   mpn_gcdext_subdiv_step, but there may be some other.

At least it's my current understanding that those changes should be
sufficient: there's nothing fundamental in the subquadratic gcdext
code that breaks canonicalization of cofactors.

I don't think I'll get the time to do this soon; it's not an awful lot
of work but I think it requires a day or so of full attention, and I
usually get hacking time only for one or a few hours at a time.

Regards,
/Niels