Is it possible to implement the extended euclidean algorithm with unsigned machine words?

Question

I'm trying to find a way to implement EEA using uint64_t in C, on a system that will not support 128-bit integers. The problem is that it seems like there is always a case where some variable or another will overflow, generating incorrect results.

I know it can be done in /signed/ machine words, and there are plenty of answers on here and elsewhere that give pseudocode for this. Can it be done with unsigned and no overflow? Or do you need to have a larger integer size available?

Answer 1

The coefficients in the Euclidian algorithm alternate in sign (after the initial zeros vanish), so we can compute them using magnitudes only and reconstructing the sign at the end from the iteration parity.

The extended Euclidian algorithm finds the multiplicative inverses a and b of positive relatively prime integers u and v with respect to each other. That is, we find a and b such that a•u is congruent to 1 modulo v and b•v is congruent to 1 modulo u. We start with two equations:

• u = 1•u + 0•v.
• v = 0•u + 1•v.

Then the algorithm is:

• Let x and y be the left sides of the latest two equations (y is later).
• If y is 1, we are done. The coefficients of u and v on the right side are the inverses of u and v, respectively.
• Otherwise, let q be the integer quotient x divided by y. Append a new equation that equals the earlier equation minus the later equation times q.

This implements the Euclidean algorithm on the left side of the equation, so we know it terminates in 1 for relatively prime u and v. Then the last equation has the form:

• 1 = a•u + b•v.

In this, it is clear a is a multiplicative inverse of u modulo v and b is a multiplicative inverse of b modulo u.

Observe that in the third equation, the coefficient of u will be 1−0•q. Thus, it will be positive. And the coefficient of v will be 0−1•q. Thus, it will be negative. In the fourth equation, the coefficient of u will be positive. From that point on, we are always subtracting a negative number from a positive number or a positive number from a negative number, and thus we alternate signs.

In the n^th equation, the coefficient of u is non-negative if n is odd and non-positive if n is even, and vice-versa for the coefficient of v.

Therefore, we can implement this with unsigned arithmetic by keeping only the magnitudes of the coefficients:

• Given magnitude g of a known-non-negative coefficient and a magnitude h of a known-non-positive coefficient, calculate the new magnitude as g+h•q.
• Given magnitude g of a known-non-positive coefficient and a magnitude h of a known-non-negative coefficient, calculate the new magnitude as g+h•q.

Example for 13 and 10:

• 13 = 1•13 + 0•10.
• 10 = 0•13 + 1•10.
• Quotient for 13/10 is 1. Compute 13−1•10 = 3, 1+0•1 = 1, 0+1•1 = 1. This gives:
• 3 = 1•13 − 1•10. (Sign is known implicitly, not computed.)
• Quotient for 10/3 is 3. Compute 10−3•3 = 1, 0+1•3 = 3, 1+1•3 = 4. This gives:
• 1 = −3•13 + 4•10.

At this point, whatever variable we are using to hold the coefficient of u (13) contains 3, but we know it is negative because we are in an even iteration (the fourth). So an inverse of u is −3. If desired, we can add u (calculated as u minus the stored magnitude) to get a positive result.

I have proven these calculations never exceed v for the coefficient of u or u for the coefficient of v but do not have access to that proof. (It is embedded in source code for a previous employer.)

Answer 2

Eric Postpischil is absolutely right, but the answer is hard to read, especially if you are looking for code that just works, so here you go:

template<typename T>
typename std::enable_if<std::is_unsigned<T>::value,tuple<T,T>>::type
constexpr extended_euclidean(const T a, const T b)
{
  T r0 = a;
  T r1 = b;
  T s0 = 1;
  T s1 = 0;
  T t0 = 0;
  T t1 = 1;
  size_t n = 0;
  while (r1) {
    T q = r0 / r1;
    r0 = r0>q*r1?r0-q*r1:q*r1-r0; swap(r0,r1);
    s0 = s0+q*s1; swap(s0,s1);
    t0 = t0+q*t1; swap(t0,t1);
    ++n;
  }
  // gcd = r0
  if (n%2) s0=b-s0;
  else     t0=a-t0;
  return tuple<T,T>({s0,t0});
}

So basically we are doing the standard euclidean algorithm, but when computing the remainders we only care about the magnitude and when updating the coefficients we can just add. The sign of the magnitude needs to be fixed in the end using the parity, using a counter n.