GCD of 2 numbers

The algorithm reduces the problem of finding the GCD by repeatedly applying these identities:

1. gcd(0, v) = v, because everything divides zero, and v is the largest number that divides v. Similarly, gcd(u, 0) = u. gcd(0, 0) is not defined.
2. If u and v are both even, then gcd(u, v) = 2·gcd(u/2, v/2), because 2 is a common divisor.
3. If u is even and v is odd, then gcd(u, v) = gcd(u/2, v), because 2 is not a common divisor. Similarly, if u is odd and v is even, then gcd(u, v) = gcd(u, v/2).
4. If u and v are both odd, and u ≥ v, then gcd(u, v) = gcd((u−v)/2, v). If both are odd and u <v, then gcd(u, v) = gcd((v-u)/2, u). These are combinations of one step of the simple Euclidean algorithm, which uses subtraction at each step, and an application of step 3 above. The division by 2 results in an integer because the difference of two odd numbers is even.
5. Repeat steps 3–4 until u = v, or (one more step) until u = 0. In either case, the result is 2^kv, where k is the number of common factors of 2 found in step 2.

Since this definition is tail-recursive, a loop can be used to replace the recursion.

unsigned int gcd(unsigned int u, unsigned int v)
{
int shift;

/* GCD(0,x) := x */
if (u == 0 || v == 0)
  return u | v;

/* Let shift := lg K, where K is the greatest power of 2
   dividing both u and v. */
for (shift = 0; ((u | v) & 1) == 0; ++shift) {
    u >>= 1;
    v >>= 1;
}

while ((u & 1) == 0)
  u >>= 1;

/* From here on, u is always odd. */
do {
    while ((v & 1) == 0)  /* Loop X */
      v >>= 1;

    /* Now u and v are both odd, so diff(u, v) is even.
       Let u = min(u, v), v = diff(u, v)/2. */
    if (u <= v) {
             v -= u;        
    } else {                      
             int diff = u - v;
             u = v; 
             v = diff;      
    } 
    v >>= 1;
    } while (v != 0);
    return u;
}