# Range for the coeffients returned by gcdext

Niels Möller nisse at lysator.liu.se
Tue Sep 16 22:26:09 CEST 2003

```Computing gcdext means that given two positive integers a, b, we
compute g = gcd(a,b) and two values u, v such that

g = u a + v b

Now, these are not uniquely determined, which can be seen from

g = u a + v b = (u + b) a + (v - a) b

In fact, u is uniquely defined mod b, and v is defined mod a.

So which representatives should be returned? It would make some sense
to always choose u such that 0 <= u < b. Then, we get 0 >= v > -a,
since

v = (g - u a) / b > - ua / b > - ba / b = - a

(except that u = 0 implies v = 1 > 0), Alternatively, one could choose

-(b-1)/2 <= u <= (b-1)/2  for odd b, and
-b/2 < u <= b/2           for even b

So what is expected from mpz_gcdext and mpn_gcdext? I find it a little
odd that mpn_gcdext is documented as returning a positive or negative
u, but most other mpn functions deal with negative numbers only when
absolutely necessary.

By the way, I think I've understood how to do binary gcdext (although
in the context of Schönhage's algorithm, I can use that only for the
one limb basecase):

Assume that A and B are odd. Through out the algorithm, we maintain two
odd numbers a and b and two values u0 and u1, 0 <= u0, u1 < B, such
that

a = u0 A + v0 B
b = u1 A + v1 B

(where we don't keep track of the v0 and v1). We initialize with

a = A, u0 = 1
b = B, u1 = 0

Then the main loop goes like

for (;;)
{
if (a == b)
{
Done, return g = a and u = u0.
}
else if (a > b)
{
swap(a, b); swap (u0, u1);
}

b -= a; /* This makes b even */
if (u1 >= u0)
u1 -= u0;
else
u1 += B - u0;

do { /* We loop until b is odd again */

/* As b = u1 A + v1 B is even, while A and B are odd,
either both or none of u1, v1 is even */

if (u1 % 2 == 0)
u1 /= 2;
else
/* u1 + B is even.
Overflow-safe way of computing (u1 + B) / 2 */
u1 = u1/2 + B/2 + 1;

b /= 2;
} while (b % 2 == 0);
}

This algorithm always returns a u >= 0.

Regards,
/Niels
```