Matt Wetmore

nn choose kk the Haskell way

Binomial coefficients are a surprisingly interesting subject of study, so interesting in fact that Sherlock Holmes’ arch-enemy Professor Moriarty began his career with a paper on them. Incidentally, this fictional Treatise on the Binomial Theorem has a longer Wikipedia page than any real work of mathematics I’ll ever write. But I digress. Today we’re going to look at a novel way to calculate binomial coefficients in Haskell, the programming language that is basically math.

The binomial coefficients (nk)n \choose k draw their name from the following “Binomial formula”:

(1+x)n=k=0n(nk)xk(1 + x)^n = \sum_{k = 0}^n {n \choose k} x^k

Here (nk)n \choose k is the coefficient on the kkth power of xx in the expansion of the left-hand side. However, a “purer” definition of (nk)n \choose k is the combinatorial one: the number of ways to choose kk objects from a set of nn. From this definition one can see why these numbers show up as the coefficients above[1]. We can also use this point of view to understand some other binomial sum identities – for example, for a finite set AA with nn elements, the cardinality of the power set P(A)\mathcal P(A) (the set of all subsets of AA) is 2n2^n. This gives a way to understand why

2n=k=0n(nk)2^n = \sum_{k=0}^n {n \choose k}

which we can of course derive from the binomial formula above but whatever. For each kk with 0kn0 \leq k \leq n, there are (nk)n \choose k ways to specify a distinct set of kk elements, i.e. a subset of size kk. So the number of subsets of any size is just the sum of these numbers over all possible kk, and there are 2n2^n total subsets[2].

Another nice identity is this one:

(nk)=nk(n1k1){n \choose k} = \frac{n}{k}{n-1 \choose k-1}

Combinatorially we can understand it like so: to count the number of ways to pick kk elements from a set of nn elements, pick one of the nn elements and look at the number of ways to pick the other k1k-1 elements from the remaining n1n-1 elements. There are nn ways to pick the first element, so that’s n(n1k1)n{n-1 \choose k-1} ways of picking subsets. But we overcounted – each subset is counted kk times, one time for each of its elements that we pick first. So after dividing by kk to fix this, we’ve shown the identity holds.

The identity is so nice because it gives a natural recursive algorithm to compute the coefficients; just keep applying the identity to get

(nk)=nkn1k1nk+11(nk0){n \choose k} = \frac{n}{k}\cdot\frac{n-1}{k-1}\cdots\frac{n-k+1}{1}{n-k\choose 0}

where (nk0)=1{n - k \choose 0} = 1. And since this is a post about (nk)n \choose k in Haskell, I’m going to turn this into a recursive Haskell function right? Unfortunately someone beat me to it, and they even used more identities than I have, so let’s take a different approach and go back to the Binomial formula. If we can expand the polynomial (1+x)n(1 + x)^n into a sum, we can just read off the coefficients for each power of xx. Easy right?

A polynomial only has finitely-many powers of xx, so we can represent a polynomial such as a0+a1x+a2x2++anxna_0 + a_1x + a_2x^2 + \ldots + a_nx^n as a list [a0,,an][a_0, \ldots, a_n]. In Haskell, we use the typeclass Num to describe things that we can add, multiply, negate, and so on; polynomials can do these things too, so let’s define a Num instance for lists of as, where a is an instance of Num, representing polynomials with coefficients in a:

instance Num a => Num [a] where
  fromInteger n = [fromInteger n]
  (x:xs) + (y:ys) = (x + y) : (xs + ys)
  xs + [] = xs
  [] + ys = ys
  (x:xs) * (y:ys) = (x*y) : ([x] * ys + xs * (y:ys))
  _ * _ = []

This is missing a few Num functions, but we won’t need them[3]. We add two polynomials by adding the respective coefficients for each power of xx, and the rather complicated expression for multiplication can be derived without much work[4].

So now [1,1][1,1] is our representation for 1+x1 + x. Exponentiation comes for free since Haskell implements it using repeated squaring. So if we try it out, [1,1]^2 == [1,2,1], and [1,1]^4 == [1,4,6,4,1]. Translating these back to math world, (1+x)2=1+2x+x2(1+x)^2 = 1 + 2x + x^2 and (1+x)4=1+4x+6x2+4x3+x4(1+x)^4 = 1 + 4x + 6x^2 + 4x^3 + x^4. So now it’s only natural to define (nk)n \choose k by:

choose :: Int -> Int -> Int
choose n k = ([1,1]^n) !! k

where xs !! k is the kkth element of xs. Now we may quickly concoct some quality coefficients:

3 `choose` 2 == 3
10 `choose` 6 == 210
434 `choose` 87 == 4614992264942942144

or just generate an infinite Pascal’s triangle with map ([1,1]^) [0..].

[[1]
,[1,1]
,[1,2,1]
,[1,3,3,1]
,[1,4,6,4,1]
,[1,5,10,10,5,1]
,[1,6,15,20,15,6,1]
,[1,7,21,35,35,21,7,1]
,[1,8,28,56,70,56,28,8,1]
,[1,9,36,84,126,126,84,36,9,1]
, ...

This is honestly an absurd[5] way to calculate binomial coefficients, but the underlying concept of manipulating polynomials as lists is very cool. We can take this in more directions, such as using infinite streams to represent power series, but perhaps that is a topic for a later post.

Acknowledgments

Thanks to Ira Kones and Rory Bokser for proofreading this post.


  1. When multiplying out (1+x)n(1+x)^n, for each of the nn clauses (1+x)(1 + x) we have (nk)n \choose k ways of choosing kk of the xx's to get xkx^k ↩︎

  2. If you’re wondering how we knew the number of subsets in the first place, note that each subset SAS \subseteq A can be identified with a function χS:A{0,1}\chi_S : A \to \{0,1\}, where χS(a)=1\chi_S(a) = 1 if and only if aSa \in S, and there are 2n2^n such functions by a simple counting argument. ↩︎

  3. You can disable warnings for this by passing ghc the flag -fno-warn-missing-methods. ↩︎

  4. Let B=(b0+b1x)B = (b_0 + b_1x\ldots). Then (a0+a1x)B(a_0 + a_1x\ldots)\cdot B =a0B+(a1x)B= \;a_0B + (a_1x\ldots)B =a0b0+a0(b1x)+a1xB= \;a_0b_0 + a_0(b_1x\ldots) + a_1xB =a0b0+x[a0(b1)+a1B]= \;a_0b_0 + x[a_0(b_1\ldots) + a_1B] ↩︎

  5. Which is why I like it so much. ↩︎