Learning Haskell and Number Theory: GCD and Higher Order Functions

This is the fourth installment of my series “Learning Haskell and Number Theory”, in which I work through Gareth and Mary Jones’ text book, Elementary Number Theory, and translate its ideas and algorithms into Haskell.  In the last installment, I wrote a function to compute the greatest common divisor of two numbers.  In this installment, we will play with higher order functions, and culminate with the extension of greatest common divisors to an arbitrary number of inputs.

The second of the two GCD implementations in the previous installment looked 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

The expression iterate next (abs a, abs b) is interesting, in that its first parameter (next) is itself a function. This general idea can be extended to all sorts of tasks, and is one of the cornerstones of functional programming.  It goes by a name: higher order functions.  (The general intuition here is that a first order function is one that operates on basic pieces of data, such as integers or strings; a second order function is one that operates on first order functions, and so on.  I don’t like this intuition because it ignores the fact that what’s chosen as “simple” data in any given language is arbitrary; but I don’t get to pick these names.)  A higher order function can take other function as parameters, and it can evaluate to a function as its result.  If it does both of these, it can be seen as a sort of transformer of functions.

I intend to write a higher-order function, and then talk about some of them that are provided in the Haskell Prelude.  Before I can get there, though, I need to clean up something that many of you pointed to in comments earlier, and that’s the meaning of the greatest common divisor of zero and zero.  The Jones text calls it undefined, which makes very literal sense.  Everything divides zero, and the set of integers has no greatest element.  Therefore, there can be no greatest divisor of zero.  However, as many of you pointed out, defining gcd(0, 0) = 0 has some nice benefits.

One of those benefits is that it allows me to continue where I want to go.  (I’m still following the book, but in fact the book is wrong here.  As an exercise, if you’ve got the book, see if you can find the error in the book’s answer to exercise 1.9.  There must be an error, but the statement for which it wants a proof is not true!)

Therefore, we’ve got:

gcd1 :: forall a. Integral a => a -> a -> a
gcd1 a b = gcd1' (abs a) (abs b)
    where gcd1' a 0 = a
          gcd1' a b = gcd1' b (a `rem` b)

Suppose we want the greatest common divisor of more than two numbers.  Following example 1.9 of the Jones text, we could define:

gcd' :: forall a. Integral a => a -> a -> a -> a
gcd' a b c = gcd (gcd a b) c

gcd'' :: forall a. Integral a => a -> a -> a -> a -> a
gcd'' a b c d = gcd (gcd (gcd a b) c) d

gcd''' :: forall a. Integral a => a -> a -> a -> a -> a -> a
gcd''' a b c d e = gcd (gcd (gcd (gcd a b) c) d) e

Oops… looks like we’re going on forever here. We could define gcd'''' or even worse.  Better is to exploit this fact:

*NumTheory> :qc \a -> gcd1 0 a == abs a
+++ OK, passed 100 tests.

Great!  So that gives us a starting point, 0, into which we can add numbers one by one until we have the greatest common divisor of them all.  We need a function that converts a list of integral values into a single integral value representing the GCD.

gcdmany :: forall a. Integral a => [a] -> a
gcdmany [] = 0
gcdmany xs = let x = init xs
                 y = last xs
             in  gcd1 (gcdmany x) y

The new functions init and last do what they say. The init function returns the initial part of a list — more precisely, all but the last element. The last function returns the last element of the list. Note that x is again a list!  There is, of course, one more thing that may be tricky about this code: it is recursive, like most useful things in Haskell.  Some careful thought, though, should convince you that this does exactly what the book suggests: from left to right, it applies the greatest common divisor algorithm to each pair of numbers until it gets to the end.

It’s quite inefficient, though.  At turns out that in haskell, it’s much faster to pick the first element from a list than the last one, because lists are singly linked in the forward dimension.  A more complicated, but faster, implementation is:

gcdmany :: forall a. Integral a => [a] -> a
gcdmany xs = gcdmany' 0 xs
    where gcdmany' p []     = p
          gcdmany' p (x:xs) = gcdmany' (gcd1 p x) xs

If you’ve got the hang of Haskell, you can check that this works.  The only thing I haven’t mentioned yet is the list constructor, :, which is basically LISP’s cons function if you know other functional languages.  If the code remains a mystery to you (which would be normal if you’re new to the language) you can study it, or just trust me.  The problem is that this is a little more complicated than we’d like to write again and again, and it turns out that you want to do exactly this lots and lots in Haskell code!  A higher order function would be handy, and the standard library provides one, called foldl, that does exactly what we’re looking for.

gcdmany :: forall a. Integral a => [a] -> a
gcdmany xs = foldl gcd1 0 xs

The last thing we can do, just be fancy, is called eta reduction.  This is the observation that if f(x) = g(x) for all x, then f = g.  This property is sometimes known as extensional equivalence, and it’s a property of pure functions, by definition. The xs at the end of the left and right-hand sides of the Haskell equation match up just like the xs in the math equation, and we can remove them.

gcdmany :: forall a. Integral a => [a] -> a
gcdmany = foldl gcd1 0

And that’s pretty simple.