Ordinary generating functions

A generating function stores a whole sequence \(a_0, a_1, a_2, \ldots\) as the coefficients of a single power series \(A(x) = \sum_{n \ge 0} a_n x^n\). Herbert Wilf put it memorably: a generating function is a clothesline on which we hang up a sequence of numbers for display. The point is that natural combinatorial operations on the objects being counted become algebraic operations on the series. Sums of disjoint classes become sums of series, independent choices become products, and recursive structure becomes an equation that can be solved. Euler used generating functions to study partitions in 1748; Laplace developed them for probability; today they are the standard language of enumeration.

Formal power series

We treat \(A(x)\) as a formal power series: an infinite list of coefficients written with powers of \(x\) as place-holders, with no claim that it converges for any particular \(x\). Formal series are added coefficient by coefficient and multiplied by the Cauchy product \[\Big(\sum a_n x^n\Big)\Big(\sum b_n x^n\Big) = \sum_n \Big(\sum_{k=0}^n a_k b_{n-k}\Big)x^n.\] Each coefficient of the product is a finite sum, so no convergence is needed. A series has a multiplicative inverse exactly when its constant term is nonzero; for example \((1 - x)(1 + x + x^2 + \cdots) = 1\), so \[\frac{1}{1 - x} = \sum_{n \ge 0} x^n.\] We write \([x^n]A(x)\) for the coefficient of \(x^n\). When a series does converge near 0, the formal and analytic answers agree, so calculus can be borrowed freely, and lesson 6 uses analysis to estimate coefficients.

A small dictionary of operations. Multiplying by \(x\) shifts a sequence right: \(xA(x)\) has coefficients \(0, a_0, a_1, \ldots\). Differentiating and multiplying by \(x\) multiplies the \(n\)th term by \(n\): \(xA'(x) = \sum na_nx^n\). Multiplying by \(1/(1 - x)\) produces partial sums: \([x^n]\frac{A(x)}{1 - x} = a_0 + \cdots + a_n\). Replacing \(x\) by \(cx\) multiplies \(a_n\) by \(c^n\).

Key expansions. From \(\frac{1}{1 - x}\) and the product rule below, \[\frac{1}{(1 - x)^k} = \sum_{n \ge 0}\binom{n + k - 1}{k - 1}x^n,\] the multiset count of lesson 2, and \(\frac{1}{1 - cx} = \sum c^nx^n\). Newton's generalized binomial theorem gives, for any number \(\alpha\), \[(1 + x)^\alpha = \sum_{n \ge 0}\binom{\alpha}{n}x^n, \qquad \binom{\alpha}{n} = \frac{\alpha(\alpha - 1)\cdots(\alpha - n + 1)}{n!}.\]

The product rule

Suppose objects of class \(\mathcal{A}\) have sizes, with \(a_n\) objects of size \(n\), and similarly \(b_n\) for class \(\mathcal{B}\). An ordered pair \((\alpha, \beta)\) with \(\alpha \in \mathcal{A}\), \(\beta \in \mathcal{B}\) has size \(|\alpha| + |\beta|\). The number of pairs of total size \(n\) is \(\sum_k a_k b_{n-k}\), which is exactly \([x^n]A(x)B(x)\). So the generating function for pairs is the product of the generating functions. This is why generating functions work: choosing several things independently and adding their sizes corresponds to multiplying series.

Making change. In how many ways can 50 cents be made from pennies, nickels, dimes and quarters? Using any number of pennies contributes \(1 + x + x^2 + \cdots = \frac{1}{1 - x}\), any number of nickels \(1 + x^5 + x^{10} + \cdots = \frac{1}{1 - x^5}\), and so on. A way to make change is an independent choice for each coin, so the answer is \[[x^{50}]\frac{1}{(1 - x)(1 - x^5)(1 - x^{10})(1 - x^{25})} = 49.\] George Pólya used this example in a 1956 article on picture-writing to show how generating functions organize a computation; the coefficient is found by multiplying truncated polynomials, as in the program below. For a dollar the count is 242, and including half-dollar coins makes it 292.

Restricted compositions. The compositions of \(n\) using parts 1 and 2 correspond to sequences of such parts, and a sequence of length \(k\) has generating function \((x + x^2)^k\). Summing over \(k\), \[F(x) = \sum_{k \ge 0}(x + x^2)^k = \frac{1}{1 - x - x^2}.\] So the counts are the Fibonacci numbers, as MATH-260 found by tiling a strip with squares and dominoes. In general, sequences of objects from a class with generating function \(P(x)\), \(P(0) = 0\), have generating function \(1/(1 - P(x))\).

Vandermonde's identity. Comparing coefficients of \(x^n\) in \((1 + x)^{m}(1 + x)^{k} = (1 + x)^{m+k}\) gives \(\sum_j \binom{m}{j}\binom{k}{n-j} = \binom{m+k}{n}\). Choosing \(n\) people from \(m\) women and \(k\) men is the combinatorial meaning.

Dice and factorization

The sum of two ordinary dice has generating function \(D(x)^2\) with \(D(x) = x + x^2 + \cdots + x^6\). Is there another pair of dice, with positive integer faces, giving the same distribution of sums? Factor \[D(x) = x(1 + x)(1 + x + x^2)(1 - x + x^2).\] Any pair of dice with the same sums has face polynomials \(P, Q\) with \(PQ = D^2\), each with \(P(1) = Q(1) = 6\) faces and no constant term. The factors at \(x = 1\) take values \(1, 2, 3, 1\), so each of \(P\) and \(Q\) needs one factor \(x\), one \(1 + x\) and one \(1 + x + x^2\), and the two factors \(1 - x + x^2\) can be split between them in only two ways. Putting both in one die gives the Sicherman dice, found by George Sicherman and published by Martin Gardner in 1978: \[P(x) = x(1 + x)(1 + x + x^2) = x + 2x^2 + 2x^3 + x^4,\] faces 1, 2, 2, 3, 3, 4, and \(Q(x) = x(1 + x)(1 + x + x^2)(1 - x + x^2)^2\), faces 1, 3, 4, 5, 6, 8. This is the only other solution.

Solving recurrences

Generating functions solve linear recurrences mechanically. Take \(a_0 = 0\), \(a_1 = 1\), \(a_n = 5a_{n-1} - 6a_{n-2}\) for \(n \ge 2\). Multiply the recurrence by \(x^n\) and sum over \(n \ge 2\): \[A(x) - x = 5x\,A(x) - 6x^2A(x),\] so \[A(x) = \frac{x}{1 - 5x + 6x^2} = \frac{x}{(1 - 2x)(1 - 3x)} = \frac{1}{1 - 3x} - \frac{1}{1 - 2x}\] by partial fractions. Hence \(a_n = 3^n - 2^n\). The denominator is the characteristic polynomial written backwards, which explains the characteristic-root method of MATH-260; repeated roots give terms like \(1/(1 - rx)^2 = \sum (n + 1)r^nx^n\), the source of the extra factor \(n\).

For the Fibonacci numbers \(F_0 = 0\), \(F_1 = 1\), the same steps give \(F(x) = \frac{x}{1 - x - x^2}\). Factor \(1 - x - x^2 = (1 - \phi x)(1 - \psi x)\) with \(\phi = \frac{1 + \sqrt5}{2}\) and \(\psi = \frac{1 - \sqrt5}{2}\). Partial fractions give \(F(x) = \frac{1}{\sqrt5}\Big(\frac{1}{1 - \phi x} - \frac{1}{1 - \psi x}\Big)\), and reading off coefficients yields Binet's formula \(F_n = (\phi^n - \psi^n)/\sqrt5\).

Inhomogeneous recurrences are no harder. For the towers of Hanoi, \(h_0 = 0\) and \(h_n = 2h_{n-1} + 1\). Summing \(h_nx^n\) over \(n \ge 1\) gives \(H(x) = 2xH(x) + \frac{x}{1 - x}\), so \(H(x) = \frac{x}{(1 - x)(1 - 2x)} = \frac{1}{1 - 2x} - \frac{1}{1 - x}\) and \(h_n = 2^n - 1\).

The Catalan generating function

Let \(C(x) = \sum C_nx^n\). The first-return decomposition of lesson 2, \(C_n = \sum_{i=0}^{n-1}C_iC_{n-1-i}\) for \(n \ge 1\), says that a nonempty Dyck path is an up-step, a path, a down-step and a path. In generating functions, \[C(x) = 1 + x\,C(x)^2.\] This quadratic has solutions \(C(x) = \frac{1 \pm \sqrt{1 - 4x}}{2x}\). Only the minus sign gives a power series, since the plus sign makes the expression blow up at \(x = 0\). So \[C(x) = \frac{1 - \sqrt{1 - 4x}}{2x}.\]

Extracting coefficients. By Newton's theorem, \([x^m]\sqrt{1 - 4x} = \binom{1/2}{m}(-4)^m\). A short computation shows \(\binom{1/2}{m}(-4)^m = -\frac{2}{m}\binom{2m - 2}{m - 1}\) for \(m \ge 1\). Then \(C_n = -\frac12[x^{n+1}]\sqrt{1 - 4x} = \frac{1}{n+1}\binom{2n}{n}\), recovering the formula of lesson 2 by pure algebra. The same method counts binary trees, triangulations and every other family with the same decomposition; this is the start of the symbolic method developed by Philippe Flajolet and Robert Sedgewick in their 2009 book Analytic Combinatorics.

from fractions import Fraction
from math import comb

def mul(a, b, N):
    out = [0] * (N + 1)
    for i, x in enumerate(a[:N + 1]):
        if x:
            for j, y in enumerate(b[:N + 1 - i]):
                out[i + j] += x * y
    return out

def geometric_in_steps(step, N):
    return [1 if n % step == 0 else 0 for n in range(N + 1)]

N = 100
ways = [1] + [0] * N
for coin in [1, 5, 10, 25]:
    ways = mul(ways, geometric_in_steps(coin, N), N)
print("ways to make 50 cents:", ways[50], "  one dollar:", ways[100])
with_half = mul(ways, geometric_in_steps(50, N), N)
print("one dollar including half-dollars:", with_half[100])

def inverse(a, N):
    inv = [Fraction(0)] * (N + 1)
    inv[0] = Fraction(1, a[0])
    for n in range(1, N + 1):
        inv[n] = -sum(a[k] * inv[n - k] for k in range(1, min(n, len(a) - 1) + 1)) / a[0]
    return inv

fib = mul([0, 1], inverse([1, -1, -1], 20), 20)
print("x/(1-x-x^2):", [int(c) for c in fib])
seq = mul([0, 1], inverse([1, -5, 6], 12), 12)
print("x/(1-5x+6x^2):", [int(c) for c in seq])
print("3^n - 2^n:     ", [3 ** n - 2 ** n for n in range(13)])

die = [0, 1, 1, 1, 1, 1, 1]
two = mul(die, die, 12)
a = [0] * 9; b = [0] * 9
for f in [1, 2, 2, 3, 3, 4]:
    a[f] += 1
for f in [1, 3, 4, 5, 6, 8]:
    b[f] += 1
print("ordinary dice sums: ", two[2:13])
print("Sicherman dice sums:", mul(a, b, 12)[2:13])

def gen_binom(alpha, m):
    r = Fraction(1)
    for i in range(m):
        r *= (alpha - i)
        r /= (i + 1)
    return r

sqrt_coeffs = [gen_binom(Fraction(1, 2), m) * (-4) ** m for m in range(13)]
catalan = [-sqrt_coeffs[n + 1] / 2 for n in range(12)]
print("Catalan from (1 - sqrt(1-4x))/(2x):", [int(c) for c in catalan])
print("matches binom(2n,n)/(n+1):", [int(c) for c in catalan] == [comb(2 * n, n) // (n + 1) for n in range(12)])
C = [int(c) for c in catalan]
lhs = C[:10]
rhs = ([1] + [0] * 9)
sq = mul(C, C, 9)
rhs = [1] + [sq[n - 1] for n in range(1, 10)]
print("C(x) = 1 + x C(x)^2 checked to x^9:", lhs == rhs)

ways to make 50 cents: 49   one dollar: 242
one dollar including half-dollars: 292
x/(1-x-x^2): [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
x/(1-5x+6x^2): [0, 1, 5, 19, 65, 211, 665, 2059, 6305, 19171, 58025, 175099, 527345]
3^n - 2^n:      [0, 1, 5, 19, 65, 211, 665, 2059, 6305, 19171, 58025, 175099, 527345]
ordinary dice sums:  [1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1]
Sicherman dice sums: [1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1]
Catalan from (1 - sqrt(1-4x))/(2x): [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786]
matches binom(2n,n)/(n+1): True
C(x) = 1 + x C(x)^2 checked to x^9: True

Multiplying truncated series answers the change problem in a fraction of a second, inverting \(1 - x - x^2\) reproduces the Fibonacci numbers, and the Sicherman dice give exactly the same distribution of sums as ordinary dice. The binomial series for \(\sqrt{1 - 4x}\) produces the Catalan numbers, and the functional equation checks term by term.