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
More information about the gmp-devel
mailing list