Range for the coeffients returned by gcdext

[Niels Möller]

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