Cofactor canonicalisation of mpn_gcdext

Niels Möller nisse at lysator.liu.se
Sun Nov 22 22:02:35 CET 2009

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

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