Cofactor canonicalisation of mpn_gcdext

[Bill Allombert]

On Fri, Nov 20, 2009 at 11:27:44AM +0G100, Torbjorn Granlund wrote:
> GMP 4.3 and earlier versions do not return the same cofactor from
> mpn_gcdext (and as a consequence also from mpz_gcdext).

First thanks a lot for tackling this issue. 

> Assume we're computing g = us + vt, where g = gcd(u,v).
> 
> Earlier GMP versions seem to have computed s such that 0 <= 2|s| <= |v|/g.
> For this to hold, s sometimes will be negative.

Yes, and the fact that s could be negative was clearly documented, so
somehow the user expected such inequality.

> GMP 4.3 computes s,t that are much larger.  It is not completely clear
> what the maximum values are.  We intend to tighten this.

This is good news!

> What criteria should we apply?
> 
> We could follow the old GMP route of 0 <= 2|s| <= |v|/g, but it might be
> better to instead promise 0 <= s <= |v|/g.  I.e., s is now non-negative.
> We might also just promise the more lax 0 <= |s| <= |v|/g.
> 
> Please express any opinions you have about this as soon as possible.

I would strive to return the exact same result than GMP 4.2 for several 
reasons:

1) backward compatibility: There are a number of software compiled with GMP
<=4.2 and due to the great ABI stability of GMP, they will automatically use
GMP 4.3 when GMP is upgraded and start to return unexpected result.
Futhermore, it seems very difficult to write code that would work both
with GMP <=4.2 and a GMP version that would return s in 0 <= s < |v|/g
because the user does not know |v|/g.

2) Allowing s == |v|/g is problematic: assumes u=a*v then g=v and we expect the
trivial solution s=0: i.e u*0+v*1=v. However the other solution s=v/g=1 is also
allowed: u*1-v*(a-1)=v.
(The bug in PARI was precisely due to the assumption that if v divides a then
s=0)

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, but a single addition by |v| can
convert s from 0 <= 2|s| <= |v|/g to 0 <= s < |v|/g.

However for all applications where the value of t is needed, then 0 <= 2|s| <=
|v|/g is better, since this preserve the symmetry between s and t, and give
smaller values for t. In this case converting s from 0 <= s <= |v|/g to 0 <=
2*|s| < |v|/g require an exact division (|v|/g) which is much more costly than
the addition.

As a suggestion a function mpn_modinv could be added that assume that g=1 and
compute the modular inverse in 0<=s<|v|
 
Cheers,
Bill.