Math from scratch, part twelve: Euclid and Bézout

Now that we have an integer library, let’s see if we can do some 2300-year-old mathematics. The Euclidean Algorithm as you probably recall dates from at least 300 BC and determines the greatest common divisor (GCD) of two natural numbers. The Extended Euclidean Algorithm solves Bézout’s Identity: given non-negative integers a and b, it finds integers x and y such that:

a * x + b * y == d

where d is the greatest common divisor of a and b. Obviously if we can write an algorithm that finds x and y then we can get d, and x and y have many useful properties themselves that we’ll explore in future episodes.

Of course if we can solve the problem for non-negative a and b then we can solve it for negative arguments as well by simply inverting the sign of the results, so let’s start with that:

public static Tuple<Integer, Integer> 
    ExtendedEuclideanDivision(this Integer a, Integer b)
{
    if (a < Integer.Zero)
    {
        var result = ExtendedEuclideanDivision(-a, b);
        return Tuple.Create(-result.Item1, result.Item2);
    }
    else if (b < Integer.Zero)
    {
        var result = ExtendedEuclideanDivision(a, -b);
        return Tuple.Create(result.Item1, -result.Item2);
    }

OK, so we know now that if we’ve made it to this point in the code, that a and b are both non-negative. We’re going to write a recursive algorithm here, so once more let’s remember what we need for a recursive algorithm:

A trivial base case that stops the recursion with a clearly correct answer.
A recursive case that reduces the problem to one or more smaller problems, solves the smaller problems, and combines their solutions to solve the larger problem.

Our recursive step will be based on always making b smaller but never negative, and therefore it will eventually get to zero. The minimal solution to Bézout’s identity when b is zero is that x is one, y is zero, and d is a. The greatest number that divides both a and zero is obviously a. So that’s our base case:

    else if (b == Integer.Zero)
        return Tuple.Create(Integer.One, Integer.Zero);

We now know that if we make it this far then 0 <= a and 1 <= b.

Now we compute the quotient q and remainder r of a and b, remembering that the definition of both is:

a = b * q + r
0 <= r < b

We need a recursive step. We recurse on b and r. Since 0 <= r < b we know that we are getting closer to or have arrived at zero, so the recursive step is in fact solving a smaller problem:

    else
    {
        var q = a / b;
        var r = a % b;
        var result = ExtendedEuclideanDivision(b, r);
        var s = result.Item1;
        var t = result.Item2;

OK, let’s stop and review for a moment. We have recursively computed s and t such that the following is true by the hypothesis that the recursive step is correct:

b * s + r * t == e

where e is the greatest common divisor of b and r. We are trying to discover x and y such that

a * x + b * y == d

where d is the greatest common divisor of a and b. My claim is that t is a suitable value for x and s - q * t is a suitable value for y. Let’s see what happens when we compute a * x + b * y given these values for x and y:

  a * x + b * y 
= a * t + b * (s - q * t)
= a * t + b * s - b * q * t
= b * s + t * (a - b * q )

But a - b * q is just r, so:

= b * s + r * t
= e

So the GCD of b and r is proposed to be the GCD of a and b! We have certainly found x and y such that a * x + b * y equals something, but how do we know that it is the GCD of a and b? Let’s continue to use d as the name of the real GCD of a and b. Does e == d work?

Let’s start by asking “does e divide a and b? We know it divides b by the hypothesis that the inductive step is correct; we got back the GCD of b and r, so clearly it divides b. We also know that a = q * b + r. Since e divides both b and r we can factor it out:

a = e * ( q * b / e + r / e )

and so e must divide a evenly too.

We now know that e is a common divisor of a and b and therefore must be less than or equal to the greatest common divisor; that’s what “greatest” means! So e <= d.

We have x and y such that a * x + b * y = e. The GCD of a and b is d, so:

a * x + b * y = e
d * (a * x / d + b * y / d) = e

and therefore d divides e. Therefore d <= e.

Put it all together and we have d <= e <= d and the only way that can be true is if e is the GCD of a and b.

Let’s not forget to return that result:

        return Tuple.Create(t, s - q * t);
    }
}

And we can now use our solution to Bézout’s Identity to find the GCD of any two integers:

public static Integer GreatestCommonDivisor(this Integer a, Integer b)
{
    var result = ExtendedEuclideanDivision(a, b);
    return a * result.Item1 + b * result.Item2;
}

And we can run a little test:

var i1 = Integer.One;
var i2 = i1 + i1;
var i3 = i2 + i1;
var i5 = i2 + i3;
var i7 = i2 + i5;
var i13 = i2 * i7 - i1;
var i182 = i2 * i7 * i13;
var i294 = i2 * i3 * i7 * i7;
Console.WriteLine(i182.GreatestCommonDivisor(i294));

and get as expected 01110, which is the GCD of 182 and 294.

Next time on FAIC: We’ll use another nice property of Bézout’s Identity to solve a real-world problem.

14 thoughts on “Math from scratch, part twelve: Euclid and Bézout”

Pingback: Dew Drop – November 5, 2013 (#1,660) | Morning Dew

The complexity of this article is exponentially related to the number of times the word “obviously” appears. 😉

Reply ↓

Eric Lippert on November 5, 2013 at 12:31 pm said:

I was once in a calculus class at UW and the professor was going through a long derivation. At one point the professor said “and obviously from x it follows that y…” and a student interrupted and said “professor, is it really obvious?”. The professor stopped, stared at the board in silence for 30 seconds, and then said “yes, it is”, and continued with the derivation.

Reply ↓
- Lumj on November 6, 2013 at 6:03 am said:
  
  Haha, that’s awesome.
  
  Reply ↓
- Douglas on November 6, 2013 at 9:33 am said:
  
  I understood, “The minimal solution to Bézout’s identity when b is zero is that x is one, y is zero, and d is a.”
  
  And I understood the code that followed.
  
  I can’t make sense of, “The greatest number that divides both a and zero is obviously a.”
  
  If you mean that a * 1 + b * 0 = a (therefore d is a), that is “obvious”.
  
  Reply ↓
  - Eric Lippert on November 6, 2013 at 10:18 am said:
    
    But we must prove not just that a * 1 + b * 0 == a but that a is the GCD of a and zero. I didn’t think this fact needed proving, but let’s prove it just for the heck of it. Suppose WOLOG a is greater than zero. We seek the largest number that divides both a and zero. There are three possibilities: the GCD of a and zero is larger than, smaller than, or equal to a. The first possibility is nonsensical; no number is divided evenly by a number larger than it. (I hope I don’t need to prove this.) That leaves two possibilities; either the greatest common divisor of a and zero is a, or it is a number smaller than a. Since a is a divisor of both a and zero, no number smaller than a can be the greatest common divisor. Therefore the second possibility is nonsensical. That leaves only the third possibility: the GCD of a and zero is a.
    
    Reply ↓
    - Douglas on November 6, 2013 at 10:48 am said:
      
      You’ll have to excuse my math illiteracy. How can there be a divisor of zero? Zero can’t be divided at all. It’s the statement of zero-divisibility that I don’t understand.
    - Eric Lippert on November 6, 2013 at 10:17 pm said:
      
      Let’s go back to definitions again. We start by once again reminding you that the definition of quotient and remainder is : for natural numbers x and y, x / y –> q, x % y –> r where x == q * y + r and 0 <= r < y. Now we define divisor: a number y is said to be a divisor of x if x % y is zero.
      
      So let’s pick two numbers, say 0 and 10. Is 10 a divisor of 0? Compute the remainder. To do so we need q and r such that 0 == q * 0 + r and 0 <= r < 10. It should be clear that if 0 == 0 * q + r then r must be 0 regardless of the value we choose for q, and therefore 10 is a divisor of 0. There is nothing special about 10; all numbers are divisors of zero because all numbers leave no remainder when zero is divided by them.
    - Joshua Lawrence on November 6, 2013 at 12:54 pm said:
      
      Douglas:
      
      Yeah, it can. If you have zero apples, you can divide your apples evenly among your a friends – everyone gets nothing.
      
      I’m using apples in honour of a maths lecturer of mine, who in a stage III multivariable calculus lecture, informed the class in the middle of a long proof that “if you have one apple, plus one apple, plus one apple, you have three apples.”
      
      Glad we got that one sorted out.
    - Douglas on November 7, 2013 at 8:02 am said:
      
      Thank you for taking the time to help me sort that out. I understand now.
      
      Does the fact that zero can be evenly divided by any number hold for irrational numbers?
    - Joshua Lawrence on November 7, 2013 at 12:30 pm said:
      
      Yeah. 0/pi = 0, which divided evenly because the right-hand side has no fractional part. However, people tend not to talk about modulus and remainder for rational or irrational numbers so setting y to an irrational number in Eric’s last comment would look a little weird.
      
      I’m sure someone somewhere has generalised the concepts of modulus and remainder to any real number, but I haven’t encountered it.
    - John Payson on November 8, 2013 at 10:36 am said:
      
      Modulus is a reasonable concept for real numbers; the most obvious example being one’s angle relative to one’s starting point if one travels a certain distance around a unit circle. The only likely sources of confusion apply with integers as well: what should it mean if either operand is negative?
      
      Perhaps the most natural way to specify a modulus function on real numbers would be to specify the value of x mod y to be:
      
      If y is positive, 0 <= result < y
      If y is negative, y < result <= 0
      If y is zero, result == 0
      
      Such a definition would probably be the most useful one to apply to floating-point numbers but for one caveat: if x is a small negative number and y is (relatively speaking) a large positive number, the mathematically-correct result of x+y may not be representable; indeed, the magnitude of the difference between y and the next smaller number may be larger than the magnitude of x. From a practical perspective, I don't think that having much more precision in the result when x is near zero than when it's near y is particularly helpful and so the lack of such precision wouldn't be a particular problem, but IEEE decided to implement fmod such that the sign of the result matches the sign of the numerator, thus preserving precision but losing the otherwise-useful relation that if (x) and (x+y) are both representable, (x mod y) and ((x+y) mod y) should be equal.