This is part three of my series, “Learning Number Theory and Haskell”, in which I work through Gareth and Mary Jones’ Elementary Number Theory and translate the ideas and concepts into Haskell. In the previous two parts, I used QuickCheck to verify some theorems and translate some ideas into Haskell code.

It is widely speculated that the first real algorithm worthy of the name was due to Euclid, when he described how to calculate the greatest common divisor of two terms. What better place to start, then in writing a first serious algorithm in Haskell? (Besides, the textbook does this next.) In the process, I’ll use “equational reasoning” (a popular word in functional programming) to improve the code. Then I’ll try something completely different to finish off

First, let’s write a type signature. Greatest common divisors can be calculated with integral values only, and maps two numbers to a result. The signature is:

```gcd1 :: forall a. Integral a => a -> a -> a
```

For a first version of the algorithm, I plan to transcribe the proof from the text as code. There are a few trivial cases involving zero gotten near the top of page 6. We can write those separately using pattern matching.

```gcd1 0 0 = undefined
gcd1 a 0 = abs a
gcd1 0 b = abs b
```

There are a few more special cases for non-zero arguments.

```gcd1 a b | a <  0    = gcd1 (-a) b
| b <  0    = gcd1 a (-b)
| a == b    = a
| a <  b    = gcd1 b a
-- Phew!
| otherwise = gcd1 b (a `rem` b)
```

Let’s test against the built-in `gcd` function from the prelude. The only subtlety here is that since `gcd 0 0` is undefined, we should avoid calling it.

```*NumTheory> :qc \a b -> (a /= 0 || b /= 0) ==> gcd1 a b == gcd a b
+++ OK, passed 100 tests.
```

We’ve got a working `gcd1`. However, it is terribly unsatisfying. By following the textbook to the letter, we ended up with a function defined with eight separate cases. Yuck!

The problem is that our number theory textbook handled a lot of stuff as special cases to get it out of the way easily for the proof. Many of those special cases aren’t all that special after all, and we can clean this up by finding equivalent cases and combining them. This is really easy to do and get right in Haskell, once you’ve got the hang of it, because Haskell functions can be read and substituted just like algebra equations. This concept is often called “equational reasoning”, and is cited as an advantage of pure functional programming. The ease of doing this is a direct consequence of not having side effects.

The last two cases cover when a < b, and a > b, respectively. However, the last case subsumes the next to last one, because if a < b (with a and b both positive), then `a `rem` b` is equivalent to `a`. That eliminates one case. Similarly, the third from last case is unnecessary. If a = b, then `gcd1 a (a `rem` b)` evaluates to `gcd1 a 0`, which in turn evaluates to `abs a`. Since we know a > 0 if the evaluation gets this far, this is the same as a. This case is also unnecessary. A similar argument proves that `gcd1 0 b = b` is unnecessary. That leaves:

```gcd1 0 0 = undefined
gcd1 a 0 = abs a
gcd1 a b | a <  0    = gcd1 (-a) b
| b <  0    = gcd1 a (-b)
| otherwise = gcd1 b (a `rem` b)
```

This is better, but it’s even shorter to promote the check for negative values of a and b and for `gcd1 0 0` to an outer function. The `where` keyword can be used to define an inner function that is available only within a particular top-level function. The result looks like this:

```gcd1 0 0 = undefined
gcd1 a b = gcd1' (abs a) (abs b)
where gcd1' a 0 = a
gcd1' a b = gcd1' b (a `rem` b)
```

Eliminating all of these special cases could be interpreted as an indication that we didn’t need them to begin with.  As an exercise, try rewriting the proof in the text without assuming anything except that a and b are non-negative and are not both zero.

Rerunning the test suite to verify:

```*NumTheory> :qc \a b -> (a /= 0 || b /= 0) ==> gcd1 a b == gcd a b
+++ OK, passed 100 tests.
```

Another implementation of the greatest common divisor is interesting (though, ultimately, it doesn’t perform as well.) One way to look at the algorithm is that we are inductively defining a list of number pairs which all have the same divisors. We start with (|a|, |b|) as the first element in the list, and then define a function to get from there to the next element. That looks like this:

```gcd2 0 0 = undefined
gcd2 a b = let next (x,y)    = (y, x `rem` y)
vals          = iterate next (abs a, abs b)
Just (ans, _) = find ((== 0) . snd) vals
in  ans
```

We first define a `next` function that takes us from one pair of numbers to the next. Then we use a Haskell standard library function, `iterate`, to build a list of values resulting from successive applications of `next` to the pair of numbers. Finally, we look for the first entry from that list that has 0 for its second number. The result is the first number of that pair. You can use QuickCheck to verify that it gives the same answer as the other two algorithms above.

This introduces an explicit data structure to represent the structure of the algorithm. Sometimes, doing this can make algorithms clearer. In this case, it’s just an interesting alternative. It involves some new concepts: an infinitely long list, a higher-order function, and another interesting use of pattern matching. If you don’t understand it yet, take some time to poke around the standard library documentation. Don’t despair, though; I’ll talk about more of these ideas in future installments.

1. Dylan Thurston / Jun 6 2007 4:16 pm

Please fix a bug in the standard prelude by deleting this case:

gcd2 0 0 = undefined

The gcd of 0 and 0 is 0. I don’t know what the Haskell gods were smoking.

2. cdsmith / Jun 6 2007 5:41 pm

The definitions for corner cases can certainly be somewhat arbitrary. However, in this case there seems to be a popular agreement that the greatest common divisor of 0 and 0 is undefined. Everything divides 0, so all integers are common divisors of 0 and 0; and the set of integers has no maximum element.

3. Cale Gibbard / Jun 6 2007 11:03 pm

I was originally going to agree with that reasoning (gcd 0 0 = undefined, because everything divides 0), but there are actually a couple good reasons to define gcd 0 0 to be 0. For one, it turns the natural numbers into a complete distributive lattice, with gcd as meet and lcm as join. Secondly, it agrees with the general definition for commutative rings: d is called a common divisor of a and b if there are elements x and y such that d*x = a and d*y = b, and if every common divisor of a and b divides d, then d is called a greatest common divisor of a and b. (In general rings, they’re not unique or guaranteed to exist, but in integral domains, they’re unique up to multiplication by units.)

4. Jeff / Jun 7 2007 10:26 am

But to quote the infallible Wikipedia: “It is useful to define gcd(0, 0) = 0 and lcm(0, 0) = 0 because then the natural numbers become a complete distributive lattice with gcd as meet and lcm as join operation. This extension of the definition is also compatible with the generalization for commutative rings.”

5. Jeff / Jun 7 2007 10:28 am

Doh, sorry for the redundant comment. That will teach me to reply to a post that’s been sitting unrefreshed on my screen all night.