6. Recurrence Relations and Generating Functions

A recurrence relation defines a sequence by saying how each term depends on earlier ones. It is the natural language both for counting problems that break into smaller instances of themselves and for the running times of recursive algorithms. This chapter shows how to set up recurrences, solve linear ones exactly with characteristic equations, bound divide-and-conquer recurrences with recursion trees and the master theorem, and use generating functions, which pack a whole sequence into a single algebraic object, to solve recurrences and to count combinations under constraints.

<a id="61-recurrence-relations"></a>

6.1 Recurrence Relations

This subchapter covers how to turn a problem into a recurrence, the two classic examples (Fibonacci and the Tower of Hanoi), and the methods that convert a recurrence into a closed formula: iteration, substitution, characteristic equations, and particular solutions.

Modeling with Recurrences

Why. Many counting questions are hard to answer head-on but easy to reduce. "How many valid strings of length n are there?" becomes "how many valid strings of length n − 1 are there, and in how many ways can each be extended?" The running time of a recursive function is, likewise, its running time on smaller inputs plus the work done in the current call. A recurrence records such a reduction exactly. Once you have one, you can compute any term, prove properties by induction (3.2), and often solve for a closed form (the rest of this chapter).

What. A recurrence relation for a sequence (a_n) is an equation that expresses a_n in terms of one or more earlier terms a_0, …, a_{n−1}, for every n ≥ n_0. The initial conditions give the terms before n_0, which the rule cannot produce. A solution is a sequence that satisfies both. By strong induction (3.2), the recurrence together with its initial conditions determines exactly one sequence. The order of a recurrence is how far back it reaches: a_n = a_{n−1} + a_{n−2} has order 2.

  recurrence:            a_n = a_{n−1} + a_{n−2}     for n ≥ 2
  initial conditions:    a_0 = 1,  a_1 = 2
  unique solution:       1, 2, 3, 5, 8, 13, 21, …

  same rule, a_0 = a_1 = 1:   1, 1, 2, 3, 5, 8, …    a different sequence

The rule alone describes a whole family of sequences. The initial conditions pick out one member of that family.

How: a recipe for setting up a recurrence.

1. Define a_n precisely in words, including what n measures. Most modeling errors come from a vague definition.
2. Split the objects of size n into cases by what happens at one end: the last symbol, the first move, the largest element. The cases must be exhaustive and mutually exclusive (sum rule, 5.1).
3. In each case, remove the end piece and recognize what remains as an object counted by an earlier term.
4. Work out enough small cases by hand to supply the initial conditions.
5. Check the first few terms against brute-force enumeration.

Example 1: codewords. A system accepts a string of decimal digits as a valid codeword if it contains an even number of 0 digits. Let E_n be the number of valid codewords of length n. Split on the last digit:

  valid string of length n
  ├── last digit is 1 … 9 (9 options): the first n − 1 digits must be VALID     → 9 · E_{n−1}
  └── last digit is 0:                  the first n − 1 digits must be INVALID   → 10^{n−1} − E_{n−1}

  E_n = 9·E_{n−1} + (10^{n−1} − E_{n−1}) = 8·E_{n−1} + 10^{n−1},      E_1 = 9

The invalid strings of length n − 1 were counted with the subtraction rule (5.1): all 10^{n−1} strings minus the valid ones.

Example 2: lines cutting the plane. Draw n lines in the plane in general position: no two parallel, no three through a common point. Let L_n be the number of regions. The n-th line crosses the n − 1 earlier lines at n − 1 distinct points, which cut it into n pieces. Each piece splits the region it passes through into two, so the new line adds exactly n regions:

  L_n = L_{n−1} + n,     L_0 = 1          1, 2, 4, 7, 11, 16, 22, 29, …

Example 3: savings with deposits. An account pays 1% interest per month, and 100 is deposited at the end of every month, starting from a balance of 0. The balance after n months is B_n = 1.01·B_{n−1} + 100, with B_0 = 0. Financial models, population models, and the cost of a loop that does a fixed amount of work per iteration all have this "multiply, then add" shape.

Common shapes and where they are solved.

Shape	Example	Solved in
a_n = c·a_{n−1} + d	Tower of Hanoi, savings account	Iteration; Linear Nonhomogeneous Recurrences
a_n = c_1 a_{n−1} + ⋯ + c_k a_{n−k}	Fibonacci, tilings	Characteristic Equations
a_n = (linear rule) + F(n)	Codewords: 8E_{n−1} + 10^{n−1}	Linear Nonhomogeneous Recurrences
T(n) = a·T(n/b) + f(n)	Merge sort	6.2
a_{n+1} = Σ a_i a_{n−i}	Catalan numbers (5.6)	Generating functions (6.3)

from itertools import product

def codewords_rec(n):
    E = 9                                        # E_1: the nine strings "1" … "9"
    for k in range(2, n + 1):
        E = 8 * E + 10 ** (k - 1)                # E_k = 8·E_{k−1} + 10^{k−1}
    return E

def codewords_brute(n):
    return sum(1 for s in product("0123456789", repeat=n) if s.count("0") % 2 == 0)

print([codewords_rec(n) for n in range(1, 6)])     # [9, 82, 756, 7048, 66384]
print([codewords_brute(n) for n in range(1, 6)])   # [9, 82, 756, 7048, 66384]

L = [1]                                          # L_0 = 1: no lines, one region
for n in range(1, 8):
    L.append(L[-1] + n)                          # the n-th line adds n regions
print(L)                                         # [1, 2, 4, 7, 11, 16, 22, 29]

Key takeaways

• A recurrence plus enough initial conditions determines a unique sequence. The rule alone does not.
• To model a problem, define a_n precisely, split by what happens at one end, and recognize smaller instances in each case.
• The cases must be exhaustive and disjoint, exactly as in the sum rule.
• Always check the first few terms of a new recurrence against brute force.

Practice

1. (Beginner) Let t_n be the number of ways to tile a 2 × n board with 1 × 2 dominoes. Find a recurrence and initial conditions for t_n, and compute t_6. Hint: look at the left edge. It is covered either by one vertical domino or by two stacked horizontal ones.
2. (Intermediate) Find a recurrence for the number of bit strings of length n that do not contain three consecutive 0s, with initial conditions. Verify it by brute force for n ≤ 12. Hint: split by the number of 0s after the last 1 (0, 1, or 2).
3. (Interview) A message of uppercase letters is encoded as A → 1, …, Z → 26 and concatenated, so "226" can be decoded as BZ, VF, or BBF. Count the decodings of a digit string of length n in O(n) time. Hint: let d_i be the number of decodings of the first i digits, and split on whether the last letter used one digit or two. Pay attention to the digit 0.

Fibonacci and Tower of Hanoi

Why. These two puzzles are the standard first examples of recurrences, and between them they show the two behaviors you meet most often. In Fibonacci, each term is the sum of the two previous terms. In the Tower of Hanoi, each term is a constant multiple of the previous term plus a constant. Both grow exponentially, but for different reasons, and both are solved exactly later in this subchapter.

Fibonacci numbers. Fibonacci's puzzle (1202) starts with one newborn pair of rabbits. A pair becomes fertile after one month and then produces one new pair every month, and rabbits never die. At month n, the pairs alive are all the pairs alive at month n − 1, plus one newborn pair for every pair that was already alive at month n − 2 (those are now fertile):

  month n:           1   2   3   4   5   6   7
  older pairs:       0   1   1   2   3   5   8      = total at month n − 1
  newborn pairs:     1   0   1   1   2   3   5      = total at month n − 2
  total f_n:         1   1   2   3   5   8  13

  f_n = f_{n−1} + f_{n−2},    f_0 = 0,  f_1 = 1     (3.3)

Fibonacci numbers appear throughout computing:

Where	Why Fibonacci
Tiling a 1 × n strip with squares and dominoes	The last tile is a square (leaving n − 1) or a domino (leaving n − 2): f_{n+1} tilings
Climbing n stairs 1 or 2 at a time	The same recurrence as the tiling problem
Bit strings of length n with no 11 (5.1)	f_{n+2} strings
Euclid's algorithm (4.2)	Consecutive Fibonacci numbers are the slowest input
AVL trees of height h (Chapter 9)	The minimum number of nodes satisfies a Fibonacci-like recurrence

The ratio f_{n+1}/f_n tends to the golden ratio φ = (1 + √5)/2 ≈ 1.618, and f_n ≈ φⁿ/√5 (proved in Characteristic Equations). So f_n has about n·log₁₀ φ ≈ 0.209·n decimal digits.

Computing f_n quickly. Plain iteration takes n additions. The matrix identity [[1, 1], [1, 0]]ⁿ = [[f_{n+1}, f_n], [f_n, f_{n−1}]] and square-and-multiply (4.3) bring this down to O(log n) arithmetic operations. Expanding the matrix squaring gives the fast doubling formulas, which avoid the matrices altogether:

  f_{2k}   = f_k · (2·f_{k+1} − f_k)
  f_{2k+1} = f_k² + f_{k+1}²

Tower of Hanoi. There are three pegs, A, B, and C, and n disks of different sizes stacked on A, largest at the bottom. The goal is to move the whole stack to C. You may move only one top disk at a time, and you may never place a disk on a smaller one.

The recursive strategy has three phases: move the top n − 1 disks out of the way to B (using C as the spare), move the largest disk to C, then move the n − 1 disks from B onto it (using A as the spare). Pegs are listed bottom to top:

  phase                    A           B           C           moves
  start                    [3 2 1]     []          []
  1. top n − 1 to B        [3]         [2 1]       []          H_{n−1}
  2. largest to C          []          [2 1]       [3]         1
  3. top n − 1 onto C      []          []          [3 2 1]     H_{n−1}

  H_n = 2·H_{n−1} + 1,    H_0 = 0          0, 1, 3, 7, 15, 31, 63, …

The strategy is optimal. Consider any solution. Before the largest disk can move for the first time, the n − 1 smaller disks must all be on the third peg, because the largest disk needs an empty target peg and nothing on top of it. Getting them there takes at least H_{n−1} moves. The largest disk moves at least once. After its last move, the n − 1 smaller disks must all end up on top of it, which takes at least another H_{n−1} moves. So every solution uses at least 2H_{n−1} + 1 moves, and the recursive strategy achieves exactly that.

The values 1, 3, 7, 15, 31 suggest H_n = 2ⁿ − 1, which the next topic proves. For the legendary 64 disks, that is about 1.8 × 10¹⁹ moves, roughly 585 billion years at one move per second.

	Fibonacci	Tower of Hanoi
Recurrence	f_n = f_{n−1} + f_{n−2}	H_n = 2H_{n−1} + 1
Order	2	1
Homogeneous	Yes	No (the + 1)
Closed form	(φⁿ − ψⁿ)/√5	2ⁿ − 1
Growth factor per step	φ ≈ 1.618	2

def hanoi(n, src="A", via="B", dst="C"):
    """Yield the moves (disk, from, to) that transfer n disks from src to dst."""
    if n == 0:
        return
    yield from hanoi(n - 1, src, dst, via)   # phase 1: park the n − 1 smaller disks on via
    yield (n, src, dst)                      # phase 2: move the largest disk
    yield from hanoi(n - 1, via, src, dst)   # phase 3: bring the n − 1 disks back on top

print(list(hanoi(2)))         # [(1, 'A', 'B'), (2, 'A', 'C'), (1, 'B', 'C')]
print([sum(1 for _ in hanoi(n)) for n in range(1, 11)])
# [1, 3, 7, 15, 31, 63, 127, 255, 511, 1023]

pegs = {"A": [3, 2, 1], "B": [], "C": []}
for disk, s, d in hanoi(3):
    assert pegs[s][-1] == disk and (not pegs[d] or pegs[d][-1] > disk)   # every move is legal
    pegs[d].append(pegs[s].pop())
print(pegs)                   # {'A': [], 'B': [], 'C': [3, 2, 1]}

def fib_pair(n):
    """Return (f_n, f_{n+1}) by fast doubling: O(log n) arithmetic operations."""
    if n == 0:
        return (0, 1)
    a, b = fib_pair(n // 2)                  # a = f_k, b = f_{k+1}, with k = n // 2
    c = a * (2 * b - a)                      # f_{2k}
    d = a * a + b * b                        # f_{2k+1}
    return (d, c + d) if n % 2 else (c, d)

print([fib_pair(n)[0] for n in range(11)])   # [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
print(len(str(fib_pair(10**4)[0])))          # 2090 digits, about 0.209 · n

Key takeaways

• f_n = f_{n−1} + f_{n−2} grows by a factor of about φ ≈ 1.618 per step. It counts tilings, stair climbs, and bit strings without 11.
• f_n can be computed in O(log n) arithmetic operations by matrix powers or fast doubling.
• Hanoi: H_n = 2H_{n−1} + 1 from the three-phase strategy, and the lower-bound argument shows that no strategy does better.
• H_n = 2ⁿ − 1: the move count doubles with every added disk.

Practice

1. (Beginner) List the 7 moves for n = 3 disks. Check that the smallest disk moves on every odd-numbered move.
2. (Intermediate) Prove that a 1 × n strip has f_{n+1} tilings by squares and dominoes. Then count the tilings with exactly k dominoes to show that f_{n+1} = Σ_k C(n − k, k), the shallow diagonals of Pascal's triangle (5.3).
3. (Interview) In a variant of Hanoi, disks may only move between adjacent pegs (A↔B and B↔C, never directly A↔C). Find the minimum number of moves for n disks from A to C. Hint: the largest disk must travel A → B → C. Count how many times the n − 1 smaller disks must move as a full stack, write a recurrence, and solve it by iteration (next topic).

Iteration and Substitution Method

Why. Before any general theory, two hands-on methods work on almost every recurrence. You can unroll the recurrence until a pattern appears (iteration), or guess the answer and confirm it by induction (substitution). They are also the tools for checking the answers of every other method in this chapter.

Iteration (backward substitution). Replace a_{n−1} with its own definition, then a_{n−2}, and so on, until you reach an initial condition. Do not simplify too early: keep the structure visible so the pattern shows. For Hanoi:

  H_n = 2H_{n−1} + 1
      = 2(2H_{n−2} + 1) + 1          = 2²H_{n−2} + 2 + 1
      = 2²(2H_{n−3} + 1) + 2 + 1     = 2³H_{n−3} + 2² + 2 + 1
      ⋮
      = 2ᵏH_{n−k} + 2^{k−1} + ⋯ + 2 + 1         after k steps
      = 2ⁿH_0 + 2^{n−1} + ⋯ + 2 + 1             take k = n
      = 2ⁿ − 1                                  H_0 = 0, geometric sum (3.1)

For the lines in the plane, iteration produces a sum directly: L_n = L_0 + 1 + 2 + ⋯ + n = 1 + n(n + 1)/2.

The same computation solves every recurrence of the form a_n = c·a_{n−1} + d:

  a_n = cⁿ·a_0 + d·(c^{n−1} + ⋯ + c + 1) = cⁿ·a_0 + d·(cⁿ − 1)/(c − 1)      if c ≠ 1
  a_n = a_0 + n·d                                                         if c = 1

This covers Hanoi (c = 2, d = 1), the savings account from the previous topic, and the arithmetic and geometric progressions of 3.1.

Forward iteration computes a_0, a_1, a_2, … and looks for a pattern: Hanoi gives 0, 1, 3, 7, 15, 31, which suggests 2ⁿ − 1. This needs a guess, but it works even when the algebra of backward iteration is messy.

Iteration's weak spot. The "⋮" line is a pattern spotted by eye, not a proof. To make the result rigorous, prove it by induction. That is exactly the substitution method.

Substitution (guess and verify).

1. Guess the form of the solution, with unknown constants.
2. Substitute the guess into the recurrence and solve for the constants that make it hold for all n.
3. Fit the remaining constants to the initial conditions.
4. Conclude by induction (3.2): step 3 is the basis, and step 2 is the inductive step.

Example. Solve T(0) = 0, T(n) = 2T(n − 1) + n. The homogeneous part doubles and the extra term is linear, so guess T(n) = A·2ⁿ + Bn + C:

  A·2ⁿ + Bn + C  =  2(A·2^{n−1} + B(n − 1) + C) + n
                 =  A·2ⁿ + (2B + 1)n + (2C − 2B)

  compare n-terms:     B = 2B + 1          ⇒  B = −1
  compare constants:   C = 2C − 2B         ⇒  C = 2B = −2
  initial condition:   T(0) = A + C = 0    ⇒  A = 2

  T(n) = 2^{n+1} − n − 2                   check: T(1) = 1, T(2) = 4, T(3) = 11

Good guesses come from small values, from iteration, from the table in Linear Nonhomogeneous Recurrences, and from recursion trees (6.2).

Pitfall: proving bounds by substitution. When you prove an inequality such as T(n) ≤ c·n, the inductive step must produce the same constant c. The following "proof" that merge sort is linear is wrong:

  Claim: T(n) = 2T(n/2) + n satisfies T(n) ≤ c·n.
  "Step": T(n) ≤ 2·c·(n/2) + n = c·n + n = (c + 1)·n = O(n)          WRONG

  The hypothesis was T(n) ≤ c·n, but the step only shows T(n) ≤ (c + 1)·n.
  The constant grows by 1 at every level of recursion, so the claim was never proved.
  (In fact T(n) = Θ(n log n), see 6.2.)

def check(recurrence, closed_form, initial, n_max=30):
    """Compare a recurrence with a conjectured closed form for n = 0 … n_max.
    recurrence(n, a) returns a_n, given the list a = [a_0, …, a_{n−1}]."""
    seq = list(initial)
    while len(seq) <= n_max:
        seq.append(recurrence(len(seq), seq))
    bad = [n for n in range(n_max + 1) if seq[n] != closed_form(n)]
    return bad[:3] or "ok"                       # the first counterexamples, or "ok"

print(check(lambda n, a: 2 * a[n - 1] + 1, lambda n: 2**n - 1, [0]))            # ok   Hanoi
print(check(lambda n, a: a[n - 1] + n, lambda n: 1 + n * (n + 1) // 2, [1]))    # ok   lines
print(check(lambda n, a: 2 * a[n - 1] + n, lambda n: 2**(n + 1) - n - 2, [0]))  # ok
print(check(lambda n, a: 2 * a[n - 1] + n, lambda n: 2**(n + 1) - n, [0]))      # [0, 1, 2]   wrong guess

A numeric check like this catches wrong guesses quickly. Agreement on 31 values is still not a proof: the induction is.

Key takeaways

• Iteration unrolls the recurrence down to the initial conditions and usually ends in a sum from 3.1.
• a_n = c·a_{n−1} + d has the solution cⁿa_0 + d(cⁿ − 1)/(c − 1) for c ≠ 1.
• Substitution: guess a form with unknown constants, match coefficients, fit initial conditions, and conclude by induction.
• When proving a bound by induction, the constant must not change between the hypothesis and the conclusion.

Practice

1. (Beginner) Solve by iteration: a_n = a_{n−1} + 3 with a_0 = 2, and a_n = 5a_{n−1} with a_0 = 2.
2. (Intermediate) Solve a_n = 3a_{n−1} + 2ⁿ, a_0 = 1 by iteration, then prove your formula by induction. (Answer: 3^{n+1} − 2^{n+1}.)
3. (Interview) A recursive function scans a whole list once and then calls itself on the list without its first element (f(lst[1:])). A colleague says it runs in linear time "because each call is one pass". Write the recurrence for its running time, solve it, and propose a linear-time rewrite. Hint: T(n) = T(n − 1) + c·n. Iterate it and recognize an arithmetic series. Note that the slice itself also costs O(n).

Linear Homogeneous Recurrences

Why. Iteration works well for first-order recurrences, but on a_n = a_{n−1} + a_{n−2} it becomes unmanageable, because each unrolling step doubles the number of terms. There is, however, a large class of recurrences with a complete, mechanical solution theory: linear homogeneous recurrences with constant coefficients. Fibonacci belongs to this class, and so do many counting problems on strings and tilings. This topic defines the class and explains the structure that makes it solvable. The next topic carries out the solution.

What. A linear homogeneous recurrence relation of degree k with constant coefficients has the form

  a_n = c_1·a_{n−1} + c_2·a_{n−2} + ⋯ + c_k·a_{n−k},       c_1, …, c_k constants,  c_k ≠ 0

Each word in the name rules something out:

• Linear: earlier terms appear only to the first power and are never multiplied together.
• Homogeneous: every term on the right contains some a_j. There is no stand-alone function of n, such as + 1 or + 2ⁿ.
• Constant coefficients: the c_i do not depend on n.
• Degree k: the recurrence reaches exactly k steps back (c_k ≠ 0). It needs k initial conditions, a_0, …, a_{k−1}, which then determine a unique solution.

Recurrence	Linear homogeneous, constant coefficients?
a_n = 3a_{n−1}	Yes, degree 1
f_n = f_{n−1} + f_{n−2}	Yes, degree 2
a_n = a_{n−5}	Yes, degree 5
a_n = a_{n−1} + a_{n−2}²	No: not linear
H_n = 2H_{n−1} + 1	No: not homogeneous (the + 1)
B_n = n·B_{n−1}	No: the coefficient n is not constant

Superposition. If (s_n) and (t_n) both satisfy the recurrence (ignoring initial conditions), then so does α·s_n + β·t_n for any constants α, β. Substituting it into the right-hand side and using linearity:

  c_1(αs_{n−1} + βt_{n−1}) + ⋯ + c_k(αs_{n−k} + βt_{n−k})
     = α(c_1 s_{n−1} + ⋯ + c_k s_{n−k}) + β(c_1 t_{n−1} + ⋯ + c_k t_{n−k})
     = α·s_n + β·t_n

So the solutions form a vector space. A solution is determined by its k initial values, so the space has dimension k. This suggests the strategy: find k independent "basic" solutions, then take the linear combination that matches the initial conditions.

Why geometric sequences are the basic solutions. Try a_n = rⁿ with r ≠ 0 in a_n = c_1 a_{n−1} + c_2 a_{n−2}:

  rⁿ = c_1·r^{n−1} + c_2·r^{n−2}      divide by r^{n−2}:      r² = c_1·r + c_2

So rⁿ is a solution exactly when r is a root of a fixed polynomial, the characteristic equation of the next topic. The reason geometric sequences work is that shifting a geometric sequence by one index just multiplies it by r, so the whole recurrence collapses to an equation for r. If you have seen linear differential equations, this is the same idea as trying e^{rt}, which differentiation multiplies by r.

Computing terms without solving. Even when the roots are ugly, the recurrence is a linear map on the "state", the vector of the last k terms. The companion matrix advances the state by one step:

  ┌ a_{n+1} ┐   ┌ c_1  c_2 ┐ ┌ a_n     ┐                 ┌ f_{n+1} ┐   ┌ 1  1 ┐ⁿ ┌ 1 ┐
  │         │ = │          │ │         │     Fibonacci:  │         │ = │      │  │   │
  └ a_n     ┘   └ 1    0   ┘ └ a_{n−1} ┘                 └ f_n     ┘   └ 1  0 ┘  └ 0 ┘

For degree k, the matrix is k × k: the first row holds c_1, …, c_k, and the rows below shift the state down by one. Raising it to the n-th power by square-and-multiply (4.3) gives a_n in O(k³ log n) arithmetic operations, which is how competitive-programming solutions handle n = 10¹⁸.

def linear_rec(coeffs, initial, n):
    """a_n for a_n = c_1·a_{n−1} + … + c_k·a_{n−k}, given initial = [a_0, …, a_{k−1}]."""
    k = len(coeffs)
    window = list(initial)                   # the last k terms, oldest first
    if n < k:
        return window[n]
    for _ in range(n - k + 1):
        nxt = sum(c * a for c, a in zip(coeffs, reversed(window)))   # c_1·a_{n−1} + … + c_k·a_{n−k}
        window = window[1:] + [nxt]
    return window[-1]

def mat_mul(X, Y, mod):
    return [[sum(X[i][t] * Y[t][j] for t in range(len(Y))) % mod for j in range(len(Y[0]))]
            for i in range(len(X))]

def linear_rec_fast(coeffs, initial, n, mod):
    """The same a_n modulo mod, in O(k³ log n) operations via the companion matrix."""
    k = len(coeffs)
    if n < k:
        return initial[n] % mod
    M = [list(coeffs)] + [[int(i == j) for j in range(k)] for i in range(k - 1)]   # companion matrix
    R = [[int(i == j) for j in range(k)] for i in range(k)]                         # identity
    e = n - (k - 1)
    while e:                                 # square-and-multiply (4.3)
        if e & 1:
            R = mat_mul(R, M, mod)
        M = mat_mul(M, M, mod)
        e >>= 1
    state = [[a] for a in reversed(initial)] # column vector (a_{k−1}, …, a_0)
    return mat_mul(R, state, mod)[0][0]

MOD = 10**9 + 7
print(linear_rec([1, 1], [0, 1], 10))                # 55   Fibonacci
print(linear_rec([1, 1, 1], [0, 0, 1], 10))          # 81   "tribonacci": each term sums the last three
assert all(linear_rec_fast([1, 1, 1], [0, 0, 1], n, MOD) == linear_rec([1, 1, 1], [0, 0, 1], n) % MOD
           for n in range(60))
print(linear_rec_fast([1, 1], [0, 1], 10**18, MOD))  # 209783453   f_{10^18} mod (10^9 + 7)

Key takeaways

• Linear homogeneous with constant coefficients: a_n = c_1 a_{n−1} + ⋯ + c_k a_{n−k}, with no extra function of n.
• A degree-k recurrence needs k initial conditions, which determine a unique solution.
• Linear combinations of solutions are solutions, so the solutions form a k-dimensional vector space.
• rⁿ is a solution exactly when r satisfies the characteristic equation.
• Without solving, the companion matrix gives a_n in O(k³ log n) operations.

Practice

1. (Beginner) Classify each recurrence, and give the degree of those that are linear homogeneous with constant coefficients: a_n = 3a_{n−1} + 4a_{n−2} + 5a_{n−3}, a_n = 2n·a_{n−1} + a_{n−2}, a_n = a_{n−1} + a_{n−4}, a_n = a_{n−2} + 1, a_n = a_{n−1}².
2. (Intermediate) Verify that 2ⁿ and (−1)ⁿ both satisfy a_n = a_{n−1} + 2a_{n−2}. Use superposition to find the solution with a_0 = 2, a_1 = 7.
3. (Interview) Count the strings of length n over {0, 1, 2} that contain no two consecutive 0s, modulo 10⁹ + 7, for n up to 10¹⁸. Explain why the O(n) dynamic program is infeasible and what replaces it. Hint: split by whether the string ends in 0 to get a degree-2 recurrence, then use its companion matrix.

Characteristic Equations

Why. The previous topic reduced linear homogeneous recurrences to finding the right geometric sequences. This topic turns that idea into a complete procedure. The result gives exact formulas such as Binet's formula for Fibonacci and, more important for algorithm analysis, the exact exponential growth rate of the sequence.

What. For a_n = c_1 a_{n−1} + ⋯ + c_k a_{n−k}, the characteristic equation is

  r^k − c_1·r^{k−1} − c_2·r^{k−2} − ⋯ − c_k = 0

and its roots are the characteristic roots.

Theorem (distinct roots). If the characteristic equation has k distinct roots r_1, …, r_k, then every solution has the form

  a_n = α_1·r_1ⁿ + α_2·r_2ⁿ + ⋯ + α_k·r_kⁿ

where the constants α_i are determined by the initial conditions.

Theorem (repeated roots). A root r of multiplicity m contributes m basic solutions, rⁿ, n·rⁿ, …, n^{m−1}·rⁿ. In general,

  a_n = Σ over distinct roots r_i of   (β_{i,0} + β_{i,1}·n + ⋯ + β_{i,m_i−1}·n^{m_i−1}) · r_iⁿ

Why the extra factor n? For a double root r of r² − c_1r − c_2, we have c_1 = 2r and c_2 = −r². Then c_1·(n − 1)r^{n−1} + c_2·(n − 2)r^{n−2} = 2(n − 1)rⁿ − (n − 2)rⁿ = n·rⁿ, so n·rⁿ is a solution too.

How: the procedure.

1. Write down the characteristic equation.
2. Find its roots and their multiplicities.
3. Write the general solution with unknown constants.
4. Use the k initial conditions to set up a k × k linear system, and solve it.

Worked example: Binet's formula.

  f_n = f_{n−1} + f_{n−2}   ⇒   r² − r − 1 = 0   ⇒   φ = (1 + √5)/2 ≈ 1.618,   ψ = (1 − √5)/2 ≈ −0.618

  general solution:   f_n = α·φⁿ + β·ψⁿ
  f_0 = 0:            α + β = 0
  f_1 = 1:            α·φ + β·ψ = 1   ⇒   α(φ − ψ) = 1   ⇒   α = 1/√5,   β = −1/√5

  Binet's formula:    f_n = (φⁿ − ψⁿ)/√5

Irrational numbers combine to give integers for every n. Since |ψ| < 1, the term ψⁿ/√5 is smaller than 1/2 in absolute value, so f_n is simply φⁿ/√5 rounded to the nearest integer.

Worked example: a repeated root. Solve a_n = 6a_{n−1} − 9a_{n−2} with a_0 = 1, a_1 = 6.

  r² − 6r + 9 = (r − 3)² = 0             double root r = 3
  general solution:   a_n = (β_0 + β_1·n)·3ⁿ
  a_0 = 1:            β_0 = 1
  a_1 = 6:            (1 + β_1)·3 = 6   ⇒   β_1 = 1

  a_n = (n + 1)·3ⁿ                       check: a_2 = 6·6 − 9·1 = 27 = 3·3²

What each kind of root contributes.

Root	Contribution to a_n	Behavior
Real r > 1	α·rⁿ	Exponential growth
r = 1	Constant α	Constant
|r| < 1	α·rⁿ	Decays to 0
Negative r	α·rⁿ	Alternating sign
Multiplicity m	(polynomial of degree m − 1)·rⁿ	Extra polynomial factor
Complex pair ρe^{±iθ}	ρⁿ(A·cos nθ + B·sin nθ)	Oscillation scaled by ρⁿ

For example, a_n = −a_{n−2} has r² + 1 = 0, so r = ±i with ρ = 1, and the sequence is periodic with period 4: a_0, a_1, −a_0, −a_1, a_0, ….

Growth rate. The root of largest absolute value, the dominant root, decides the growth. If r_1 is the unique dominant root and its coefficient is nonzero, then a_n ~ α_1·r_1ⁿ (the ratio tends to 1). For example, bit strings with no 11 are counted by Fibonacci numbers, so each extra bit multiplies their number by about φ ≈ 1.618, not 2. Such a string therefore carries only about log₂ φ ≈ 0.694 bits of information per symbol, which is the kind of calculation used to design constrained codes for storage media.

Floating point. Binet's formula is exact in real arithmetic, but a 64-bit float stores only about 16 significant digits, so evaluating it in floating point fails once f_n has more digits than that:

from math import sqrt

phi = (1 + sqrt(5)) / 2
psi = (1 - sqrt(5)) / 2

def fib_exact(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

def fib_binet(n):
    return round((phi**n - psi**n) / sqrt(5))    # exact in real arithmetic only

print(next(n for n in range(200) if fib_binet(n) != fib_exact(n)))   # 71: the first wrong value

def solve_order2(c1, c2, a0, a1):
    """Closed form of a_n = c1·a_{n−1} + c2·a_{n−2}, for real characteristic roots."""
    disc = c1 * c1 + 4 * c2                  # discriminant of r² − c1·r − c2
    if disc == 0:
        r = c1 / 2                           # double root
        b0, b1 = a0, a1 / r - a0             # a_n = (b0 + b1·n)·rⁿ
        return lambda n: (b0 + b1 * n) * r**n
    r1, r2 = (c1 + sqrt(disc)) / 2, (c1 - sqrt(disc)) / 2
    alpha = (a1 - a0 * r2) / (r1 - r2)       # from α + β = a0 and α·r1 + β·r2 = a1
    beta = a0 - alpha
    return lambda n: alpha * r1**n + beta * r2**n

f = solve_order2(6, -9, 1, 6)
print([round(f(n)) for n in range(6)])       # [1, 6, 27, 108, 405, 1458] = (n + 1)·3ⁿ
g = solve_order2(1, 1, 0, 1)
print([round(g(n)) for n in range(8)])       # [0, 1, 1, 2, 3, 5, 8, 13]

Use closed forms to understand growth, and integer methods (iteration, fast doubling, companion matrices) to compute exact values.

Key takeaways

• Characteristic equation: r^k − c_1r^{k−1} − ⋯ − c_k = 0. Each root r gives a basic solution rⁿ.
• A root of multiplicity m gives rⁿ, n·rⁿ, …, n^{m−1}rⁿ.
• Procedure: find the roots, write the general solution, then solve a linear system from the initial conditions.
• Binet's formula: f_n = (φⁿ − ψⁿ)/√5, so f_n is φⁿ/√5 rounded.
• The dominant root gives the exponential growth rate.

Practice

1. (Beginner) Solve a_n = 5a_{n−1} − 6a_{n−2} with a_0 = 1, a_1 = 0.
2. (Intermediate) Solve a_n = 2a_{n−1} + a_{n−2} − 2a_{n−3} with a_0 = 3, a_1 = 6, a_2 = 0. Hint: r = 1 is a root of the characteristic polynomial, so factor out (r − 1).
3. (Interview) A colleague computes the number of stair-climbing sequences for n = 10⁶ with Binet's formula in double precision and gets a meaningless number. Later they also need the result modulo 10⁹ + 7. Explain both problems and propose an exact method that is fast enough. Hint: how many significant digits does a double hold, and how many digits does f_{10⁶} have? Then recall fast doubling and the companion matrix.

Linear Nonhomogeneous Recurrences

Why. Many recurrences that come up in practice have an extra term that does not involve the sequence. Hanoi has + 1, the codewords have + 10^{n−1}, and a recursive function that does linear work per call has + n. Such a term makes the recurrence nonhomogeneous, so the characteristic-root method no longer applies directly. It still works after one extra step: find a single particular solution, then add the general solution of the homogeneous part.

What. A linear nonhomogeneous recurrence with constant coefficients has the form

  a_n = c_1·a_{n−1} + ⋯ + c_k·a_{n−k} + F(n),         F(n) not identically 0

Dropping F(n) gives the associated homogeneous recurrence.

Theorem. If a_n^{(p)} is any one solution of the nonhomogeneous recurrence (a particular solution), then every solution has the form

  a_n = a_n^{(p)} + a_n^{(h)},       where a_n^{(h)} solves the associated homogeneous recurrence

Proof. If a_n and a_n^{(p)} both solve the nonhomogeneous recurrence, subtract the two equations. The F(n) terms cancel, so a_n − a_n^{(p)} solves the homogeneous recurrence. ∎

Analogy. This is the same structure as a linear system Ax = b: all solutions are one particular solution plus anything in the null space of A.

Finding a particular solution: undetermined coefficients. When F(n) is a polynomial times an exponential, the particular solution has the same shape:

F(n)	Trial a_n^{(p)} if s is not a characteristic root	If s is a root of multiplicity m
Constant d (s = 1)	C	C·n^m
Polynomial of degree t (s = 1)	P_t(n)	n^m·P_t(n)
d·sⁿ	C·sⁿ	C·n^m·sⁿ
(polynomial of degree t)·sⁿ	P_t(n)·sⁿ	n^m·P_t(n)·sⁿ

Here P_t(n) stands for a polynomial of degree t with unknown coefficients. The extra factor n^m handles resonance: when F already looks like a homogeneous solution, the plain trial would cancel out completely.

How: the procedure.

1. Solve the associated homogeneous recurrence, keeping the constants unknown.
2. Choose a trial particular solution from the table, substitute it, and solve for its coefficients.
3. Add the two: a_n = a_n^{(p)} + a_n^{(h)}.
4. Only now fit the initial conditions, to the sum. Fitting them to the homogeneous part alone is the most common mistake.

Example 1: Hanoi again. H_n = 2H_{n−1} + 1, H_0 = 0.

  homogeneous:   r = 2                ⇒   H_n^{(h)} = α·2ⁿ
  particular:    F = 1 = 1ⁿ, and 1 is not a root   ⇒   try C:   C = 2C + 1   ⇒   C = −1
  sum:           H_n = α·2ⁿ − 1
  H_0 = 0:       α − 1 = 0   ⇒   α = 1          H_n = 2ⁿ − 1

Example 2: codewords. E_n = 8E_{n−1} + 10^{n−1}, E_1 = 9.

  homogeneous:   α·8ⁿ
  particular:    F = (1/10)·10ⁿ, and 10 is not a root   ⇒   try C·10ⁿ:
                 C·10ⁿ = 8C·10^{n−1} + 10^{n−1}   ⇒   10C = 8C + 1   ⇒   C = 1/2
  sum:           E_n = α·8ⁿ + 10ⁿ/2
  E_1 = 9:       8α + 5 = 9   ⇒   α = 1/2       E_n = (8ⁿ + 10ⁿ)/2        check: E_2 = (64 + 100)/2 = 82

Example 3: resonance. a_n = 2a_{n−1} + 2ⁿ, a_0 = 0. The plain trial C·2ⁿ fails: it gives C·2ⁿ = C·2ⁿ + 2ⁿ, that is, 0 = 2ⁿ. Since s = 2 is a root of multiplicity 1, try C·n·2ⁿ instead:

  C·n·2ⁿ = 2·C(n − 1)·2^{n−1} + 2ⁿ   ⇒   Cn = C(n − 1) + 1   ⇒   C = 1
  a_n = α·2ⁿ + n·2ⁿ,   a_0 = 0 ⇒ α = 0       a_n = n·2ⁿ          check: a_2 = 2·2 + 4 = 8

Sums are nonhomogeneous recurrences. The partial sum S_n = 1² + ⋯ + n² satisfies S_n = S_{n−1} + n². Its characteristic root is 1 and F(n) = n²·1ⁿ resonates with it, so the trial is n·(An² + Bn + C), a cubic. Matching coefficients reproduces n(n + 1)(2n + 1)/6 (3.1).

Example 4: savings. B_n = 1.01·B_{n−1} + 100, B_0 = 0. The homogeneous part is α·1.01ⁿ. The constant trial gives C = 1.01C + 100, so C = −10 000. From B_0 = 0, α = 10 000, so B_n = 10 000·(1.01ⁿ − 1).

def iterate(step, a0, n_max):
    """[a_0, …, a_{n_max}] for the first-order rule a_n = step(n, a_{n−1})."""
    seq = [a0]
    for n in range(1, n_max + 1):
        seq.append(step(n, seq[-1]))
    return seq

# Codewords, with E_0 = 1 (the empty string has zero 0s, an even number)
E = iterate(lambda n, prev: 8 * prev + 10 ** (n - 1), 1, 20)
assert E == [(8**n + 10**n) // 2 for n in range(21)]

# Resonance: a_n = 2a_{n−1} + 2ⁿ
a = iterate(lambda n, prev: 2 * prev + 2**n, 0, 20)
assert a == [n * 2**n for n in range(21)]

# Sum of squares as a nonhomogeneous recurrence
S = iterate(lambda n, prev: prev + n * n, 0, 20)
assert S == [n * (n + 1) * (2 * n + 1) // 6 for n in range(21)]

# Savings: 100 per month at 1% monthly interest, for 10 years
B = iterate(lambda n, prev: 1.01 * prev + 100, 0, 120)
print(round(B[120], 2), round(10_000 * (1.01**120 - 1), 2))   # 23003.87 23003.87

Key takeaways

• General solution = one particular solution + the general solution of the associated homogeneous recurrence.
• For F(n) = (polynomial)·sⁿ, try the same shape, multiplied by n^m if s is a characteristic root of multiplicity m.
• Fit the initial conditions only after adding the particular and homogeneous parts.
• Summations such as Σ i² are nonhomogeneous recurrences with characteristic root 1.

Practice

1. (Beginner) Solve a_n = 2a_{n−1} + 3, a_0 = 1.
2. (Intermediate) For the recurrence a_n = 6a_{n−1} − 9a_{n−2} + F(n), write the form of the trial particular solution when F(n) is n², 2ⁿ, 3ⁿ, and n·3ⁿ. Then fully solve the case F(n) = 2ⁿ with a_0 = a_1 = 0.
3. (Interview) A loan of L at monthly interest rate r is repaid with equal payments P at the end of every month, so the balance satisfies B_n = (1 + r)B_{n−1} − P with B_0 = L. Implement monthly_payment(L, r, N) in O(1) time so that the loan is cleared after exactly N payments. What changes when r = 0? Hint: the particular solution is a constant. Set B_N = 0 and solve for P. For r = 0, the characteristic root is 1, which resonates with the constant term. <a id="62-divide-and-conquer-recurrences"></a>

6.2 Divide-and-Conquer Recurrences

Algorithms that split a problem into a few pieces of a fraction of the original size lead to recurrences in which n is divided rather than decremented. This subchapter introduces such algorithms, then two ways to solve their recurrences: drawing the recursion tree, and applying the master theorem.

Divide-and-Conquer Algorithms

Why. If each step of an algorithm removes only one element, it needs n levels of work. If each step instead cuts the input into parts that are a constant fraction of the size, it needs only log n levels. This is how merge sort beats the quadratic sorts and how binary search beats a linear scan. The running times of such algorithms satisfy recurrences of a new shape, T(n) = a·T(n/b) + f(n), which the methods of 6.1 do not handle directly.

What. A divide-and-conquer algorithm has three steps:

1. Divide the input of size n into a subproblems, each of size n/b.
2. Conquer each subproblem recursively. Inputs below a fixed size are solved directly.
3. Combine the subresults into the answer.

If the divide and combine steps take f(n) time together, the running time satisfies

  T(n) = a · T(n/b) + f(n)
         │     │       └── work outside the recursive calls (divide + combine)
         │     └────────── each subproblem is a 1/b fraction of the input
         └──────────────── number of recursive calls

We use Θ(g(n)) to mean "grows at the same rate as g(n), up to constant factors". Chapter 10 defines it formally (10.2).

Algorithm	Recurrence	Solution
Binary search	T(n) = T(n/2) + O(1)	Θ(log n)
Fast exponentiation (4.3)	T(n) = T(n/2) + O(1)	Θ(log n) multiplications
Merge sort	T(n) = 2T(n/2) + Θ(n)	Θ(n log n)
Finding both min and max	T(n) = 2T(n/2) + 2	3n/2 − 2 comparisons
Closest pair of points	T(n) = 2T(n/2) + Θ(n)	Θ(n log n)
Integer multiplication, 4 half-size products	T(n) = 4T(n/2) + Θ(n)	Θ(n²), no gain
Karatsuba multiplication	T(n) = 3T(n/2) + Θ(n)	Θ(n^{log₂ 3}) ≈ Θ(n^{1.585})
Strassen matrix multiplication	T(n) = 7T(n/2) + Θ(n²)	Θ(n^{log₂ 7}) ≈ Θ(n^{2.807})

Floors and ceilings. Strictly, merge sort splits n into ⌈n/2⌉ and ⌊n/2⌋. For asymptotic analysis it is safe to assume that n is a power of b: for the non-decreasing running times that appear in practice, the answer for all other n lies between the answers for the neighboring powers of b.

Solving exactly when n is a power of b. The substitution n = bᵏ turns a divide-and-conquer recurrence into an ordinary one in k, which 6.1 can solve.

  Binary search:  T(n) = T(n/2) + 1,  T(1) = 1.        Let t_k = T(2ᵏ):
                  t_k = t_{k−1} + 1,  t_0 = 1    ⇒   t_k = k + 1    ⇒   T(n) = log₂ n + 1

  Merge sort:     T(n) = 2T(n/2) + n,  T(1) = 0.       Let t_k = T(2ᵏ):
                  t_k = 2t_{k−1} + 2ᵏ,  t_0 = 0  ⇒   t_k = k·2ᵏ     ⇒   T(n) = n·log₂ n

The merge sort case is exactly the resonance example from Linear Nonhomogeneous Recurrences: the extra term 2ᵏ has the same base as the characteristic root 2, and that is where the factor log n comes from.

Karatsuba's trick. To multiply two n-digit numbers, split each into halves: x = x_1·10^m + x_0 and y = y_1·10^m + y_0. The schoolbook expansion needs four half-size products: x_1y_1, x_1y_0, x_0y_1, and x_0y_0. Karatsuba observed that the middle term can be computed from a single extra product:

  x·y = z_2·10^{2m} + z_1·10^m + z_0,  where
        z_2 = x_1·y_1
        z_0 = x_0·y_0
        z_1 = x_1·y_0 + x_0·y_1 = (x_1 + x_0)(y_1 + y_0) − z_2 − z_0     ← one multiplication, not two

Three recursive multiplications instead of four change the recurrence from 4T(n/2) + Θ(n), which is Θ(n²), to 3T(n/2) + Θ(n), which is Θ(n^{1.585}). Python's own int multiplication uses Karatsuba for large numbers.

import random
from math import log2

def merge_sort(a, stats):
    if len(a) <= 1:
        return a
    mid = len(a) // 2
    left, right = merge_sort(a[:mid], stats), merge_sort(a[mid:], stats)   # conquer: 2 calls, size n/2
    out, i, j = [], 0, 0
    while i < len(left) and j < len(right):                                # combine: Θ(n) merge
        if left[i] <= right[j]:
            out.append(left[i]); i += 1
        else:
            out.append(right[j]); j += 1
    out += left[i:] + right[j:]
    stats["moved"] += len(out)               # elements written by this merge
    return out

random.seed(1)
data = [random.randrange(10**6) for _ in range(1024)]
stats = {"moved": 0}
assert merge_sort(data, stats) == sorted(data)
print(stats["moved"], 1024 * 10)             # 10240 10240: exactly n·log₂ n for n = 2^10

def karatsuba(x, y, stats):
    """Multiply non-negative integers with 3 recursive products instead of 4."""
    if x < 10 or y < 10:
        stats["base"] += 1                   # base case: a single-digit factor
        return x * y
    m = max(len(str(x)), len(str(y))) // 2
    x1, x0 = divmod(x, 10**m)                # x = x1·10^m + x0
    y1, y0 = divmod(y, 10**m)
    z2 = karatsuba(x1, y1, stats)
    z0 = karatsuba(x0, y0, stats)
    z1 = karatsuba(x1 + x0, y1 + y0, stats) - z2 - z0    # = x1·y0 + x0·y1
    return z2 * 10**(2 * m) + z1 * 10**m + z0

for digits in (64, 128, 256, 512):
    x = random.randrange(10**(digits - 1), 10**digits)
    y = random.randrange(10**(digits - 1), 10**digits)
    stats = {"base": 0}
    assert karatsuba(x, y, stats) == x * y
    print(digits, stats["base"], digits**2)  # base-case products vs. schoolbook digit products
# 64 1045 4096
# 128 3207 16384
# 256 9723 65536
# 512 28655 262144       doubling the digits multiplies the work by about 3, not 4

Key takeaways

• Divide-and-conquer: a subproblems of size n/b, plus f(n) work to divide and combine, gives T(n) = a·T(n/b) + f(n).
• Assume n is a power of b. Substituting n = bᵏ turns the recurrence into a 6.1-style recurrence in k.
• Saving one subproblem (Karatsuba: 3 instead of 4, Strassen: 7 instead of 8) lowers the exponent of the running time.

Practice

1. (Beginner) Ternary search on a sorted array compares the target with the elements at one third and two thirds of the range, then recurses into one of the three parts. Write its recurrence for the number of comparisons and solve it exactly for n = 3ᵏ.
2. (Intermediate) To find both the minimum and the maximum of n numbers, split the list in half, solve both halves, and combine with 2 comparisons. Solve T(2) = 1, T(n) = 2T(n/2) + 2 for n = 2ᵏ, and compare the result with the 2n − 3 comparisons of the obvious approach.
3. (Interview) Count the inversions of an array (pairs i < j with A[i] > A[j]) in O(n log n) time. Write the recurrence for your algorithm. Hint: modify merge sort. When an element of the right half is output before the remaining elements of the left half, it forms an inversion with each of them.

Recursion Trees

Why. Solving T(n) = a·T(n/b) + f(n) exactly is often unnecessary, and with uneven splits it can be hard. A recursion tree is a picture of all the work the algorithm does. It shows at a glance whether most of the cost is at the root (the top-level combine), spread evenly over the levels, or at the leaves (the many tiny base cases). That tells you the answer, or at least a guess to prove by substitution (6.1).

What. In a recursion tree, each node stands for one call and is labeled with the work that call does outside its recursive calls. Its children are the calls it makes. The total running time is the sum of all labels, which is easiest to compute level by level:

  level          nodes          size of each     cost of the level
  0              1              n                f(n)
  1              a              n/b              a · f(n/b)
  2              a²             n/b²             a² · f(n/b²)
  ⋮
  i              aⁱ             n/bⁱ             aⁱ · f(n/bⁱ)
  ⋮
  log_b n        n^{log_b a}    1                Θ(n^{log_b a})        (leaves)

The tree has depth log_b n, and its number of leaves is a^{log_b n} = n^{log_b a}. Both sides equal b^{(log_b a)(log_b n)}. Remember this exponent log_b a: it reappears in the master theorem.

Example: T(n) = 2T(n/2) + n (merge sort). Every level costs exactly n, and there are log₂ n + 1 levels:

  0                           n                         n
                      ┌───────┴───────┐
  1                  n/2             n/2                2 · n/2 = n
                  ┌───┴───┐       ┌───┴───┐
  2              n/4     n/4     n/4     n/4            4 · n/4 = n
                  ⋮       ⋮       ⋮       ⋮
  log₂ n     1   1   1   1   ⋯   1   1   1   1          n · 1 = n   (n leaves)
                                                        ──────────────────────
                                  total:  n · (log₂ n + 1) = Θ(n log n)

The three shapes. The level costs often form a geometric series, and the ratio of that series decides which part of the tree dominates:

Level costs	Example	Dominant part	Total
Decreasing geometric	2T(n/2) + n²: n², n²/2, n²/4, …	Root	Θ(f(n)) = Θ(n²)
Constant	2T(n/2) + n: n, n, n, …	All levels equally	Θ(f(n) · log n) = Θ(n log n)
Increasing geometric	3T(n/2) + n: n, 3n/2, 9n/4, …	Leaves	Θ(n^{log₂ 3})

In the increasing case, the sum is dominated by its last term:

  Σ_{i=0}^{log₂ n − 1} (3/2)ⁱ · n   =   n · ((3/2)^{log₂ n} − 1)/(3/2 − 1)   =   2n · (n^{log₂ 3 − 1} − 1)   =   2n^{log₂ 3} − 2n
  plus n^{log₂ 3} leaves of cost 1:     T(n) = 3n^{log₂ 3} − 2n = Θ(n^{log₂ 3})

Uneven splits. Recursion trees also handle recurrences with subproblems of different sizes, which the master theorem cannot. Take T(n) = T(n/3) + T(2n/3) + n:

                          n                         n
               ┌──────────┴──────────┐
              n/3                  2n/3             n/3 + 2n/3 = n
          ┌────┴────┐          ┌─────┴─────┐
         n/9      2n/9       2n/9        4n/9       n
          ⋮                                ⋮        ≤ n once the short branches end

Each full level costs n. The shortest root-to-leaf path (always taking n/3) has log₃ n levels, and the longest (always taking 2n/3) has log_{3/2} n levels. So the total is between n·log₃ n and n·log_{3/2} n, which means Θ(n log n). The same argument shows that quicksort runs in Θ(n log n) even if every pivot splits the input only 1 : 9.

From picture to proof. A recursion tree is a convincing argument, but it skips details such as floors and the exact number of levels. The rigorous route is to read a guess off the tree, then prove it by substitution (6.1), keeping the constant fixed in the inductive step.

def level_costs(a, b, f, n):
    """Cost of each level of the recursion tree of T(n) = a·T(n/b) + f(n), down to size 1."""
    costs, nodes, size = [], 1, n
    while size >= 1:
        costs.append(nodes * f(size))        # aⁱ nodes, each doing f(n/bⁱ) work
        nodes, size = nodes * a, size / b
    return costs

for label, a, b, f in [("2T(n/2) + n^2", 2, 2, lambda s: s * s),
                       ("2T(n/2) + n",   2, 2, lambda s: s),
                       ("3T(n/2) + n",   3, 2, lambda s: s)]:
    c = level_costs(a, b, f, 2**10)
    print(f"{label:14} levels {[round(x) for x in c[:4]]} ... {round(c[-1])}   total {round(sum(c))}")
# 2T(n/2) + n^2  levels [1048576, 524288, 262144, 131072] ... 1024   total 2096128   ≈ 2 · root
# 2T(n/2) + n    levels [1024, 1024, 1024, 1024] ... 1024   total 11264             = n · 11 levels
# 3T(n/2) + n    levels [1024, 1536, 2304, 3456] ... 59049   total 175099           ≈ 3 · leaves

Key takeaways

• Level i of the tree for a·T(n/b) + f(n) costs aⁱ·f(n/bⁱ). There are log_b n levels and n^{log_b a} leaves.
• Decreasing level costs: the root dominates. Equal level costs: multiply by log n. Increasing level costs: the leaves dominate.
• Uneven splits: bound the total by (cost per level) × (length of the longest path).
• Use the tree to guess, and substitution to prove.

Practice

1. (Beginner) Draw the first three levels of the recursion tree for T(n) = 4T(n/2) + n. Find the cost of level i, the number of leaves, and the order of growth of T(n).
2. (Intermediate) Use recursion trees to find the order of growth of T(n) = T(n/4) + T(3n/4) + n and of T(n) = T(n/4) + T(n/2) + n. Why do the answers differ? Hint: compare the sum of the subproblem sizes with n in each case.
3. (Interview) Solve T(n) = 2T(√n) + log₂ n. Hint: substitute m = log₂ n and define S(m) = T(2^m). The new recurrence for S is one you have already seen in this subchapter.

Master Theorem

Why. Drawing a tree for every recurrence is slow, and most divide-and-conquer recurrences have the same shape. The master theorem packages the three outcomes of the recursion tree into a rule you can apply in seconds. It is the standard tool for analyzing divide-and-conquer algorithms, and it is expected knowledge in algorithm interviews.

Intuition. The master theorem is a tug-of-war between the root and the leaves. The root costs f(n). The leaves together cost Θ(n^{log_b a}), so log_b a is called the critical exponent. Whichever is polynomially larger wins. If they are the same size, every level contributes equally, and a factor log n appears.

What: simplified version. For f(n) = Θ(n^d):

  T(n) = a·T(n/b) + Θ(n^d),      a ≥ 1,  b > 1,  d ≥ 0

            ┌  Θ(n^d)              if  a < b^d     (d > log_b a: the root dominates)
  T(n)  =   ┤  Θ(n^d · log n)      if  a = b^d     (d = log_b a: all levels equal)
            └  Θ(n^{log_b a})      if  a > b^d     (d < log_b a: the leaves dominate)

Comparing a with b^d avoids computing logarithms.

Why it works. With f(n) = n^d, level i costs aⁱ·(n/bⁱ)^d = n^d·(a/b^d)ⁱ. The level costs form a geometric series with ratio q = a/b^d:

• q < 1: the series is dominated by its first term, n^d.
• q = 1: each of the log_b n + 1 levels costs n^d, for a total of Θ(n^d log n).
• q > 1: the series is dominated by its last term, n^d·q^{log_b n} = n^{log_b a}.

What: general version. Let a ≥ 1, b > 1, and f(n) ≥ 0. Compare f(n) with n^{log_b a}:

Case	Condition on f(n)	T(n)
1	f(n) = O(n^{log_b a − ε}) for some ε > 0	Θ(n^{log_b a})
2	f(n) = Θ(n^{log_b a} · logᵏ n) for some k ≥ 0	Θ(n^{log_b a} · log^{k+1} n)
3	f(n) = Ω(n^{log_b a + ε}) for some ε > 0, and a·f(n/b) ≤ c·f(n) for some c < 1 and large n	Θ(f(n))

Here O means "at most a constant times" and Ω means "at least a constant times" (10.2). The condition a·f(n/b) ≤ c·f(n) in case 3, called the regularity condition, says that level costs really do shrink geometrically. It holds automatically for polynomial f.

Applications.

Recurrence	a, b, d	Compare a with b^d	Result
Binary search: T(n/2) + 1	1, 2, 0	1 = 2⁰	Θ(log n)
Tree traversal: 2T(n/2) + 1	2, 2, 0	2 > 2⁰	Θ(n)
Merge sort: 2T(n/2) + n	2, 2, 1	2 = 2¹	Θ(n log n)
Karatsuba: 3T(n/2) + n	3, 2, 1	3 > 2¹	Θ(n^{log₂ 3})
2T(n/2) + n²	2, 2, 2	2 < 2²	Θ(n²)
Strassen: 7T(n/2) + n²	7, 2, 2	7 > 2²	Θ(n^{log₂ 7})
Naive block matrix product: 8T(n/2) + n²	8, 2, 2	8 > 2²	Θ(n³)
2T(n/2) + n log n	general version, case 2, k = 1		Θ(n log² n)

When the master theorem does not apply.

Recurrence	Problem	Use instead
T(n) = T(n − 1) + n	Subtracts from n instead of dividing it	6.1: Θ(n²)
T(n) = T(n/3) + T(2n/3) + n	Subproblems of different sizes	Recursion tree, or the Akra–Bazzi theorem
T(n) = 2T(n/2) + n/log n	f is smaller than n, but not polynomially smaller	Recursion tree: Θ(n log log n)
T(n) = n·T(n/2) + 1	a is not a constant	Iteration

from math import log, log2

def power(e):
    """Format n^e, writing n^0 as 1 and n^1 as n."""
    if abs(e - round(e)) < 1e-9:
        return {0: "1", 1: "n"}.get(round(e), f"n^{round(e)}")
    return f"n^{e:.3f}"

def master(a, b, d):
    """Order of growth of T(n) = a·T(n/b) + Θ(n^d)."""
    crit = log(a, b)                         # critical exponent log_b a
    if abs(d - crit) < 1e-9:                 # a = b^d, compared with a float tolerance
        return "Θ(log n)" if d == 0 else f"Θ({power(d)} log n)"
    return f"Θ({power(max(d, crit))})"       # the larger exponent wins

for name, a, b, d in [("binary search", 1, 2, 0), ("tree traversal", 2, 2, 0),
                      ("merge sort", 2, 2, 1), ("Karatsuba", 3, 2, 1),
                      ("2T(n/2) + n^2", 2, 2, 2), ("Strassen", 7, 2, 2)]:
    print(f"{name:15} {master(a, b, d)}")
# binary search   Θ(log n)
# tree traversal  Θ(n)
# merge sort      Θ(n log n)
# Karatsuba       Θ(n^1.585)
# 2T(n/2) + n^2   Θ(n^2)
# Strassen        Θ(n^2.807)

def T(n, a, b, d):                           # exact values for n a power of b, with T(1) = 1
    return 1 if n <= 1 else a * T(n // b, a, b, d) + n**d

for k in (5, 10, 20, 30):
    n = 2**k
    print(k, round(T(n, 3, 2, 1) / n**log2(3), 4))   # the ratio tends to the hidden constant 3
# 5 2.7366
# 10 2.9653
# 20 2.9994
# 30 3.0

The last loop is a reminder of what Θ hides. The master theorem says that T(n) grows like n^{log₂ 3}. The exact solution, 3n^{log₂ 3} − 2n from the recursion-tree topic, also gives the constant 3.

Key takeaways

• Compare f(n) with n^{log_b a}, the total cost of the leaves. The larger one determines T(n). If they are equal, multiply by log n.
• Simplified version for f(n) = Θ(n^d): compare a with b^d.
• Case 3 needs the regularity condition, which polynomial f satisfies automatically.
• The theorem does not cover subtractive recurrences, uneven splits, or f that falls in the gap between cases. Use iteration or recursion trees for those.

Practice

1. (Beginner) Apply the master theorem to T(n) = 4T(n/2) + n, 4T(n/2) + n², 4T(n/2) + n³, and 9T(n/3) + n.
2. (Intermediate) Explain why the master theorem does not apply to T(n) = 2T(n/2) + n/log₂ n. Then use a recursion tree to show that T(n) = Θ(n log log n). Hint: level i costs n/(log₂ n − i), and the sum is related to the harmonic numbers (3.1).
3. (Interview) You can design a divide-and-conquer algorithm in three ways. A: 5 subproblems of half the size, with linear combine time. B: 2 subproblems of half the size, with quadratic combine time. C: 9 subproblems of a third of the size, with quadratic combine time. Which is asymptotically fastest? Hint: for each design, compare a with b^d and write the resulting exponent. A logarithmic factor breaks a tie between equal exponents. <a id="63-generating-functions"></a>

6.3 Generating Functions

A generating function encodes an entire sequence as the coefficients of a single power series, so that operations on sequences become algebra on functions. This subchapter introduces ordinary and exponential generating functions and uses them to solve recurrences, including some the characteristic equation cannot handle, and to count combinations under constraints.

Ordinary Generating Functions

Why. A sequence is infinitely many numbers. Working with it term by term is what makes recurrences hard. A generating function packs the whole sequence into one object. Herbert Wilf described it as "a clothesline on which we hang up a sequence of numbers for display". The payoff is that shifting a sequence, summing it, or convolving it with another becomes multiplication or division of functions. Recurrences then turn into algebraic equations, and counting problems turn into polynomial multiplication.

What. The ordinary generating function (OGF) of a sequence a_0, a_1, a_2, … is the power series

  G(x) = a_0 + a_1·x + a_2·x² + ⋯ = Σ_{n≥0} a_n·xⁿ

The notation [xⁿ] G(x) means "the coefficient of xⁿ in G(x)", so [xⁿ] G(x) = a_n. A finite sequence has a polynomial as its generating function. For example, the binomial theorem (5.3) says that (1 + x)^m is the generating function of C(m, 0), C(m, 1), …, C(m, m).

Formal power series. In combinatorics, x is only a placeholder that marks the position of each term, and we never substitute a number for it. Generating functions are therefore treated as formal power series: two series are equal when all their coefficients are equal, and they are added and multiplied like polynomials. Every coefficient of a product depends on only finitely many terms, so no question of convergence ever arises. In this setting, 1/(1 − x) means the unique series B(x) with (1 − x)·B(x) = 1:

  (1 − x)(1 + x + x² + x³ + ⋯)  =    1 + x + x² + x³ + ⋯
                                   − x − x² − x³ − ⋯
                                =    1                     ⇒   1/(1 − x) = Σ xⁿ

In general, a series A(x) has a multiplicative inverse exactly when its constant term a_0 is nonzero.

A table of basic generating functions.

a_n	First terms	G(x)
1	1, 1, 1, 1, …	1/(1 − x)
cⁿ	1, c, c², c³, …	1/(1 − cx)
1 for n ≤ m, then 0	1, 1, …, 1, 0, 0, …	(1 − x^{m+1})/(1 − x)
C(m, n)	1, m, C(m, 2), …	(1 + x)^m
n + 1	1, 2, 3, 4, …	1/(1 − x)²
n	0, 1, 2, 3, …	x/(1 − x)²
C(n + k − 1, n)	1, k, C(k + 1, 2), …	1/(1 − x)^k
1 if k divides n, else 0	1, 0, …, 0, 1, 0, …	1/(1 − x^k)
1/n for n ≥ 1	0, 1, 1/2, 1/3, …	ln(1/(1 − x))

The extended binomial theorem. For any real number u, define C(u, k) = u(u − 1)⋯(u − k + 1)/k!, with C(u, 0) = 1. Then

  (1 + x)^u = Σ_{k≥0} C(u, k)·xᵏ

For a non-negative integer u, the series stops after x^u, and this is the ordinary binomial theorem. For u = −n, C(−n, k) = (−1)ᵏ·C(n + k − 1, k), so

  1/(1 − x)ⁿ = Σ_{k≥0} C(n + k − 1, k)·xᵏ

which is stars and bars (5.2) in one line. With u = 1/2, the theorem gives the series of √(1 − 4x), which leads to the Catalan numbers later in this subchapter.

from fractions import Fraction
from math import comb

N = 10                                       # keep the first N coefficients of every series

def series(*coeffs):
    return [Fraction(c) for c in coeffs] + [Fraction(0)] * (N - len(coeffs))

def mul(A, B):
    """Product of two truncated series: coefficient n is Σ A[i]·B[n − i]."""
    return [sum(A[i] * B[n - i] for i in range(n + 1)) for n in range(N)]

def inv(A):
    """Coefficients of 1/A(x). Requires A[0] != 0."""
    B = [1 / A[0]]
    for n in range(1, N):
        B.append(-sum(A[i] * B[n - i] for i in range(1, n + 1)) / A[0])   # forces (A·B)_n = 0
    return B

geo = inv(series(1, -1))                             # 1/(1 − x)
print([int(c) for c in geo])                         # [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
print([int(c) for c in mul(geo, geo)])               # [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]   1/(1 − x)²
print([int(c) for c in inv(series(1, -3))])          # [1, 3, 9, 27, 81, 243, 729, 2187, 6561, 19683]
cube = mul(mul(geo, geo), geo)
print(all(cube[n] == comb(n + 2, 2) for n in range(N)))   # True: 1/(1 − x)³ ↔ C(n + 2, 2)
print([int(c) for c in inv(series(1, -1, -1))])      # [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]   1/(1 − x − x²)

The last line is a preview: 1/(1 − x − x²) generates the Fibonacci numbers (shifted by one). The denominator is the Fibonacci recurrence in disguise, as the topic Solving Recurrences with Generating Functions explains.

Key takeaways

• The OGF of (a_n) is Σ a_n xⁿ, and [xⁿ] extracts a coefficient.
• Formal power series are manipulated algebraically. Convergence plays no role.
• 1/(1 − cx) = Σ cⁿxⁿ is the most important entry in the table. Most other entries follow from it.
• Extended binomial theorem: (1 + x)^u = Σ C(u, k)xᵏ for any real u. In particular, 1/(1 − x)ⁿ generates C(n + k − 1, k).

Practice

1. (Beginner) Find closed forms for the generating functions of 2, 2, 2, …, of 0, 0, 1, 1, 1, …, of 1, −1, 1, −1, …, and of 1, 0, 1, 0, 1, ….
2. (Intermediate) Find the sequences generated by 1/(1 − x)³ and by x²/(1 − 2x). Then use the extended binomial theorem to show that [xᵏ] (1 − 4x)^{−1/2} = C(2k, k).
3. (Interview) Implement a function that returns the first N coefficients of 1/A(x), given the coefficients of A(x). What is its running time, and why must a_0 be nonzero? Hint: the coefficients b_n of the inverse satisfy Σ_{i=0}^{n} a_i b_{n−i} = 0 for every n ≥ 1, which you can solve for b_n. Look up Newton iteration for power series for a faster method.

Operations on Generating Functions

Why. Generating functions are useful because natural operations on sequences correspond to simple operations on functions. Once you know these correspondences, you can build the generating function of a complicated sequence from simple ones, and read off closed forms that would otherwise take a clever induction proof.

What. Let A(x) = Σ a_n xⁿ and B(x) = Σ b_n xⁿ.

Operation on sequences	Resulting sequence	Generating function
Linear combination	α·a_n + β·b_n	α·A(x) + β·B(x)
Shift right by k (prepend k zeros)	0, …, 0, a_0, a_1, …	xᵏ·A(x)
Shift left by 1 (drop a_0)	a_1, a_2, a_3, …	(A(x) − a_0)/x
Scaling by cⁿ	cⁿ·a_n	A(cx)
Multiplying by n	n·a_n	x·A′(x)
Differences	a_n − a_{n−1}	(1 − x)·A(x)
Partial sums	a_0 + a_1 + ⋯ + a_n	A(x)/(1 − x)
Convolution	Σ_{i=0}^{n} a_i·b_{n−i}	A(x)·B(x)

Convolution is the key entry. When two series are multiplied, the coefficient of xⁿ collects one product a_i·b_j for every way of writing n = i + j:

  [x³] A(x)·B(x) = a_0·b_3 + a_1·b_2 + a_2·b_1 + a_3·b_0
                   └─ one term for every split 3 = i + j ─┘

Partial sums are a special case: they are the convolution with 1, 1, 1, …, whose generating function is 1/(1 − x). Convolution is also why generating functions count. If a_i counts the ways to do something with i units and b_j counts the ways to do something else with j units, then Σ a_i b_{n−i} counts the ways to do both with n units in total. This is the subject of Counting with Generating Functions.

Example: the sum of squares. Start from 1/(1 − x) and apply the table:

  1/(1 − x)                       = Σ xⁿ                         a_n = 1
  x · d/dx:   x/(1 − x)²          = Σ n·xⁿ                       multiply by n
  x · d/dx:   x(1 + x)/(1 − x)³   = Σ n²·xⁿ                      multiply by n again
  ÷ (1 − x):  x(1 + x)/(1 − x)⁴   = Σ (0² + 1² + ⋯ + n²)·xⁿ      partial sums

  [xⁿ] x/(1 − x)⁴  = C(n + 2, 3)          [xⁿ] x²/(1 − x)⁴ = C(n + 1, 3)      (1/(1 − x)⁴ ↔ C(n + 3, 3), then shift)

  0² + 1² + ⋯ + n² = C(n + 2, 3) + C(n + 1, 3) = n(n + 1)(2n + 1)/6

This derives the formula of 3.1 without having to guess it first.

N = 12

def mul(A, B):        return [sum(A[i] * B[n - i] for i in range(n + 1)) for n in range(N)]
def x_deriv(A):       return [n * a for n, a in enumerate(A)]           # x·A′(x)  ↔  n·a_n
def partial_sums(A):  return mul(A, [1] * N)                            # A(x)/(1 − x)
def shift(A, k):      return [0] * k + A[:N - k]                        # xᵏ·A(x)
def scale(A, c):      return [c**n * a for n, a in enumerate(A)]        # A(cx)  ↔  cⁿ·a_n

ones = [1] * N                                                          # 1/(1 − x)
squares = x_deriv(x_deriv(ones))                                        # Σ n²·xⁿ
print(squares)                  # [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121]
print(partial_sums(squares))    # [0, 1, 5, 14, 30, 55, 91, 140, 204, 285, 385, 506]
assert partial_sums(squares) == [n * (n + 1) * (2 * n + 1) // 6 for n in range(N)]
print(scale(ones, -1))          # [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1]   1/(1 + x)
print(shift(ones, 2))           # [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]         x²/(1 − x)

Key takeaways

• Shifts multiply or divide by powers of x, scaling by cⁿ substitutes cx for x, and multiplying by n applies x·d/dx.
• Partial sums divide by (1 − x), and differences multiply by (1 − x).
• Multiplying generating functions convolves the sequences. That is the bridge to counting.
• Combining a few table entries derives closed forms such as Σ i² mechanically.

Practice

1. (Beginner) Use the table to find the generating function of a_n = 3n + 2.
2. (Intermediate) Find the generating functions of a_n = n·2ⁿ and of the harmonic numbers H_n = 1 + 1/2 + ⋯ + 1/n (with H_0 = 0).
3. (Interview) The probability distribution of the sum of two fair dice is the convolution of two uniform distributions. Explain this using generating functions, and explain how to compute the distribution of the sum of 1000 dice efficiently. Hint: the generating function of one die is (x + x² + ⋯ + x⁶)/6. Raise it to a power by repeated squaring, with an FFT for each polynomial product.

Solving Recurrences with Generating Functions

Why. Characteristic equations solve linear recurrences with constant coefficients and nothing else. Generating functions solve those too, with the initial conditions built in from the start, and they also handle recurrences that involve convolutions, such as the Catalan recurrence, which no characteristic equation can touch. The method is always the same: turn the recurrence into an equation for G(x), solve it, and expand.

How: the method.

1. Let G(x) = Σ a_n xⁿ.
2. Multiply both sides of the recurrence by xⁿ and sum over every n for which the recurrence holds.
3. Express each sum in terms of G(x) using the shift rules. The initial conditions supply the missing low-order terms.
4. Solve the resulting equation for G(x).
5. Split G(x) into partial fractions and read off the coefficients, using 1/(1 − rx) = Σ rⁿxⁿ.

Example 1: a first-order recurrence. Solve a_n = 3a_{n−1} + 2 for n ≥ 1, with a_0 = 1.

  Σ_{n≥1} a_n·xⁿ   =  3·Σ_{n≥1} a_{n−1}·xⁿ  +  2·Σ_{n≥1} xⁿ
  G(x) − 1         =  3x·G(x)               +  2x/(1 − x)

  G(x)·(1 − 3x) = 1 + 2x/(1 − x) = (1 + x)/(1 − x)
  G(x) = (1 + x)/((1 − x)(1 − 3x)) = A/(1 − x) + B/(1 − 3x)

  1 + x = A(1 − 3x) + B(1 − x):     x = 1:    2 = −2A       ⇒   A = −1
                                    x = 1/3:  4/3 = 2B/3    ⇒   B = 2

  G(x) = 2/(1 − 3x) − 1/(1 − x)     ⇒     a_n = 2·3ⁿ − 1

Example 2: Fibonacci. The recurrence f_n = f_{n−1} + f_{n−2} holds for n ≥ 2, with f_0 = 0 and f_1 = 1.

  Σ_{n≥2} f_n·xⁿ = Σ_{n≥2} f_{n−1}·xⁿ + Σ_{n≥2} f_{n−2}·xⁿ
  F(x) − x       = x·F(x)             + x²·F(x)

  F(x) = x/(1 − x − x²) = x/((1 − φx)(1 − ψx)) = (1/√5)·(1/(1 − φx) − 1/(1 − ψx))

  ⇒   f_n = (φⁿ − ψⁿ)/√5          (Binet's formula again, with φ + ψ = 1 and φψ = −1)

The denominator 1 − x − x² is the characteristic polynomial r² − r − 1 with its coefficients reversed, and its roots are the reciprocals 1/φ and 1/ψ of the characteristic roots. This is true in general: a degree-k linear recurrence with constant coefficients has a rational generating function P(x)/Q(x) with Q(x) = 1 − c_1x − ⋯ − c_k xᵏ. The two methods are the same theory seen from two sides. This connection also gives a quick asymptotic rule: the root of Q closest to 0 determines the exponential growth rate.

Example 3: Catalan numbers. Here generating functions do something characteristic equations cannot. The Catalan recurrence (5.6) is a convolution:

  C_{n+1} = Σ_{i=0}^{n} C_i·C_{n−i},   C_0 = 1

  Σ_{n≥0} C_{n+1}·x^{n+1} = x · Σ_{n≥0} (Σ_i C_i·C_{n−i})·xⁿ
  C(x) − 1                = x · C(x)²                               convolution ↔ product

  x·C(x)² − C(x) + 1 = 0   ⇒   C(x) = (1 − √(1 − 4x)) / (2x)

The quadratic formula gives two roots. With + in front of the square root, the expression would contain a 1/x term and would not be a power series, so the minus sign is the right one. Expanding √(1 − 4x) = Σ C(1/2, k)·(−4x)ᵏ with the extended binomial theorem and simplifying gives C_n = C(2n, n)/(n + 1).

	Characteristic equation	Generating functions
Handles	Linear recurrences with constant coefficients	Also convolutions, systems of recurrences, some variable coefficients
Initial conditions	Fitted at the end	Built in from the start
Main work	Finding roots, solving a linear system	Algebra and partial fractions
Result	A closed form	A closed form, plus a generating function that can be reused for further counting and asymptotics

from fractions import Fraction
from math import comb

def coeffs_of_ratio(P, Q, n_terms):
    """First n_terms coefficients of P(x)/Q(x), by power-series division (requires Q[0] = 1)."""
    out = []
    for n in range(n_terms):
        c = (P[n] if n < len(P) else 0) - sum(Q[i] * out[n - i] for i in range(1, min(n, len(Q) - 1) + 1))
        out.append(c)                        # from (Q·G)_n = P_n
    return out

print(coeffs_of_ratio([1, 1], [1, -4, 3], 8))     # (1 + x)/((1 − x)(1 − 3x))
# [1, 5, 17, 53, 161, 485, 1457, 4373]  =  2·3ⁿ − 1
print(coeffs_of_ratio([0, 1], [1, -1, -1], 10))   # x/(1 − x − x²)
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

def gen_binom(u, k):
    """Extended binomial coefficient C(u, k) for any rational u."""
    out = Fraction(1)
    for i in range(k):
        out *= (u - i) / Fraction(i + 1)
    return out

sqrt_coeffs = [gen_binom(Fraction(1, 2), k) * (-4) ** k for k in range(12)]   # √(1 − 4x)
catalan = [-sqrt_coeffs[n + 1] / 2 for n in range(11)]   # [xⁿ] (1 − √(1 − 4x))/(2x)
print([int(c) for c in catalan])     # [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796]
assert catalan == [comb(2 * n, n) // (n + 1) for n in range(11)]

Key takeaways

• Method: multiply by xⁿ, sum, rewrite in terms of G(x), solve, use partial fractions, and read off coefficients.
• Linear recurrences with constant coefficients have rational generating functions, whose denominators are the reversed characteristic polynomials.
• Convolution recurrences, such as Catalan's, become polynomial equations for G(x).
• The root of the denominator closest to 0 gives the exponential growth rate.

Practice

1. (Beginner) Use generating functions to solve a_n = 2a_{n−1} for n ≥ 1, with a_0 = 3.
2. (Intermediate) Use generating functions to solve a_n = 5a_{n−1} − 6a_{n−2} with a_0 = 1, a_1 = 0, and check that you get the same answer as with the characteristic equation.
3. (Interview) Show that the generating function for the number of bit strings of length n with no two consecutive 0s is (1 + x)/(1 − x − x²). Without solving for the coefficients exactly, what is their exponential growth rate? Hint: set up the recurrence by the last bit and apply the method. For the growth rate, find the root of the denominator closest to 0 and take its reciprocal.

Counting with Generating Functions

Why. Many counting problems ask how many ways there are to reach a total: distribute 8 cookies among 3 children with each getting 2 to 4, make change for a dollar, roll a sum of 10 with three dice. Stars and bars (5.2) handles only the simplest constraints. Generating functions handle arbitrary constraints on each part with the same method: write one factor per part and multiply.

What. The exponent of x records how many units are used, so x^e stands for "this part uses e units". For each part, write a polynomial or series whose terms are the allowed amounts, then multiply:

  number of solutions of  e_1 + e_2 + ⋯ + e_k = n,  e_i ∈ S_i      =      [xⁿ]  Π_i ( Σ_{s ∈ S_i} xˢ )

Expanding the product picks one term from each factor in every possible way, and each pick contributes x^{e_1 + ⋯ + e_k}. So the coefficient of xⁿ counts exactly the picks whose exponents add up to n.

  8 cookies, 3 children, each child gets 2, 3, or 4:

     child 1           child 2           child 3
  (x² + x³ + x⁴)  ·  (x² + x³ + x⁴)  ·  (x² + x³ + x⁴)

  picking x³ · x² · x³ = x⁸  ↔  the distribution (3, 2, 3)

  [x⁸](x² + x³ + x⁴)³ = [x²](1 + x + x²)³ = 6

Factors for common constraints.

Constraint on one part	Factor
Any amount e ≥ 0	1/(1 − x)
At least m	x^m/(1 − x)
At most m	(1 − x^{m+1})/(1 − x)
An even amount	1/(1 − x²)
A multiple of k (coins of value k)	1/(1 − xᵏ)
Used at most once	1 + x
An item of size k, used at most once	1 + xᵏ

Examples.

• Stars and bars: with no constraints, [xⁿ] 1/(1 − x)ᵏ = C(n + k − 1, n) (5.2).
• Making change: the number of ways to make 100 cents from pennies, nickels, dimes, and quarters is [x¹⁰⁰] 1/((1 − x)(1 − x⁵)(1 − x¹⁰)(1 − x²⁵)) = 242. The coin-change dynamic program of 5.6 is this product, computed one factor at a time.
• Dice: the number of ways to roll a sum of s with three dice is [xˢ](x + x² + ⋯ + x⁶)³. For s = 10 it is 27, so the probability is 27/216.
• Integer partitions (5.6): a partition of n chooses how many parts of each size k to use, so Σ p(n)xⁿ = Π_{k≥1} 1/(1 − xᵏ). Partitions into distinct parts use each size at most once: Π_{k≥1} (1 + xᵏ).

Euler's theorem in two lines. The theorem from 5.6 says that partitions into distinct parts are as numerous as partitions into odd parts. Using 1 + xᵏ = (1 − x^{2k})/(1 − xᵏ):

  Π_{k≥1} (1 + xᵏ) = Π_{k≥1} (1 − x^{2k})/(1 − xᵏ) = Π_{k odd} 1/(1 − xᵏ)

Every numerator 1 − x^{2k} cancels the denominator with even exponent 2k, and only the odd denominators remain. The right-hand side is the generating function for partitions into odd parts.

def poly_mul(A, B, limit):
    """Product of two polynomials, truncated after x^limit."""
    C = [0] * (limit + 1)
    for i, a in enumerate(A[:limit + 1]):
        if a:
            for j, b in enumerate(B[:limit + 1 - i]):
                C[i + j] += a * b
    return C

def factor(exponents, limit):
    """The polynomial Σ x^e over the allowed amounts e (up to limit)."""
    P = [0] * (limit + 1)
    for e in exponents:
        if e <= limit:
            P[e] += 1
    return P

def count(constraints, n):
    """Number of solutions of e_1 + … + e_k = n with each e_i in constraints[i]."""
    G = [1] + [0] * n
    for allowed in constraints:
        G = poly_mul(G, factor(allowed, n), n)   # one factor per part
    return G[n]

print(count([range(2, 5)] * 3, 8))                              # 6     cookies: 3 children, 2 … 4 each
print(count([range(0, 101, c) for c in (1, 5, 10, 25)], 100))   # 242   change for a dollar
print(count([range(1, 7)] * 3, 10), 6**3)                       # 27 216   P(sum = 10) = 27/216

n = 30                                                          # Euler: distinct parts vs. odd parts
distinct = count([(0, k) for k in range(1, n + 1)], n)          # part size k used 0 times or once
odd = count([range(0, n + 1, k) for k in range(1, n + 1, 2)], n)
print(distinct, odd)                                            # 296 296

Key takeaways

• One factor per part, and one term x^e per allowed amount e. The answer is the coefficient of xⁿ in the product.
• Products of generating functions count independent choices whose sizes add up.
• Coin change and partitions use factors 1/(1 − xᵏ). "At most once" uses 1 + xᵏ.
• Identities such as Euler's theorem become short manipulations of generating functions.

Practice

1. (Beginner) In how many ways can you choose 12 pieces of fruit from apples, bananas, and oranges, if you take at least 2 apples, at most 3 bananas, and an even number of oranges? Write the generating function and extract the coefficient.
2. (Intermediate) Show with generating functions that the number of ways to make n cents from 1-cent and 2-cent coins is ⌊n/2⌋ + 1. Then write the generating function for coins of 1, 2, and 5 cents, and compute the number of ways to make 20 cents.
3. (Interview) Find two six-sided dice, labeled with positive integers and different from the standard pair, whose sum has exactly the same distribution as the sum of two standard dice (Sicherman dice). Hint: factor x + x² + ⋯ + x⁶ = x(1 + x)(1 + x + x²)(1 − x + x²) and redistribute the factors of its square between two dice. Each die must have 6 faces (its polynomial equals 6 at x = 1) and no face labeled 0.

Exponential Generating Functions

Why. The product of two ordinary generating functions splits a total n into amounts i + j, but it does not track which items go where. That is right for identical items such as coins and cookies, but wrong for labeled objects, such as the positions in a string, the people in a seating arrangement, or the elements of a set partition. Splitting n labeled items into a group of size k and a group of size n − k can be done in C(n, k) ways, and the ordinary product misses that factor. Exponential generating functions build it in.

What. The exponential generating function (EGF) of (a_n) is

  Â(x) = Σ_{n≥0} a_n · xⁿ/n!

The factorials make products behave like a binomial convolution:

  Â(x)·B̂(x) = Σ_n ( Σ_k a_k/k! · b_{n−k}/(n − k)! ) xⁿ = Σ_n ( Σ_k C(n, k)·a_k·b_{n−k} ) · xⁿ/n!

In words: to build a labeled structure of size n from two parts, choose which k labels go to the first part (C(n, k) ways), build the first part on them (a_k ways), and build the second part on the rest (b_{n−k} ways).

a_n	EGF
1	eˣ
cⁿ	e^{cx}
n!	1/(1 − x)
1 if n is even, else 0	cosh x = (eˣ + e^{−x})/2
1 if n is odd, else 0	sinh x = (eˣ − e^{−x})/2
1 if n ≥ 1, else 0	eˣ − 1
Onto functions from an n-set to a k-set: k!·S(n, k) (5.6)	(eˣ − 1)ᵏ
Derangements D_n (5.5)	e^{−x}/(1 − x)
Bell numbers B_n (5.6)	e^{eˣ − 1}

Example 1: strings with letter-count constraints. A string of length n over an alphabet is determined by deciding which positions hold each letter, so it is a labeled structure. Each letter contributes one factor Σ x^j/j! over its allowed counts j. For example, count the DNA strings of length n (over A, C, G, T) that contain an even number of As:

  A: even count     → cosh x = (eˣ + e^{−x})/2
  C, G, T: any      → eˣ · eˣ · eˣ = e^{3x}

  Â(x) = (eˣ + e^{−x})/2 · e^{3x} = (e^{4x} + e^{2x})/2     ⇒     a_n = (4ⁿ + 2ⁿ)/2

The answer checks out for n = 1: there are (4 + 2)/2 = 3 such strings, namely C, G, and T.

Example 2: derangements. Every permutation of n elements splits uniquely into its set of fixed points and a derangement of the remaining elements. So n! = Σ_k C(n, k)·1·D_{n−k}, which is a binomial convolution. In EGFs:

  1/(1 − x) = eˣ · D̂(x)    ⇒    D̂(x) = e^{−x}/(1 − x)    ⇒    D_n = n!·Σ_{i=0}^{n} (−1)ⁱ/i!

The last step reads off the coefficients as partial sums (dividing by 1 − x) of the series of e^{−x}. This is the inclusion–exclusion result of 5.5, obtained without inclusion–exclusion.

Example 3: set partitions. A partition into exactly k blocks is a set of k non-empty blocks. Each block has EGF eˣ − 1, so labeled sequences of k blocks have EGF (eˣ − 1)ᵏ. Dividing by k! forgets the order of the blocks, which gives S(n, k). Summing over all k:

  Σ_k (eˣ − 1)ᵏ/k! = e^{eˣ − 1}       the EGF of the Bell numbers

	Ordinary (OGF)	Exponential (EGF)
Form	Σ a_n·xⁿ	Σ a_n·xⁿ/n!
Product means	Σ a_k·b_{n−k}: split an unlabeled total	Σ C(n, k)·a_k·b_{n−k}: split a labeled set
Typical objects	Multisets, coin change, integer partitions	Strings by letter counts, permutations, set partitions
"Any amount" factor	1/(1 − x)	eˣ

from fractions import Fraction
from math import factorial

N = 10

def egf(seq):  return [Fraction(a, factorial(n)) for n, a in enumerate(seq)]    # a_n → a_n/n!
def mul(A, B): return [sum(A[i] * B[n - i] for i in range(n + 1)) for n in range(N)]
def counts(A): return [int(c * factorial(n)) for n, c in enumerate(A)]          # back to a_n

exp_cx = lambda c: egf([c**n for n in range(N)])          # e^{cx}
even = egf([1 - n % 2 for n in range(N)])                 # cosh x: an even number of copies

print(counts(mul(even, exp_cx(3))))      # [1, 3, 10, 36, 136, 528, 2080, 8256, 32896, 131328]
print([(4**n + 2**n) // 2 for n in range(N)])             # the same list

derangements = mul(exp_cx(-1), egf([factorial(n) for n in range(N)]))    # e^{−x} · 1/(1 − x)
print(counts(derangements))              # [1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496]

e_minus_1 = egf([0] + [1] * (N - 1))                      # one non-empty block: eˣ − 1
bell, power = [Fraction(0)] * N, egf([1] + [0] * (N - 1)) # power holds (eˣ − 1)ᵏ, starting at k = 0
for k in range(N):
    bell = [b + p / factorial(k) for b, p in zip(bell, power)]   # add (eˣ − 1)ᵏ / k!
    power = mul(power, e_minus_1)
print(counts(bell))                      # [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147]

Key takeaways

• EGF: Σ a_n xⁿ/n!. Use it when the objects are labeled and splitting them means choosing which labels go where.
• Product of EGFs ↔ binomial convolution Σ C(n, k)a_k b_{n−k}.
• For strings with letter-count constraints, use one factor per letter: eˣ for any count, cosh x for even, sinh x for odd, eˣ − 1 for at least one.
• Derangements (e^{−x}/(1 − x)), onto functions ((eˣ − 1)ᵏ), and Bell numbers (e^{eˣ − 1}) all follow from the same product rule.

Practice

1. (Beginner) Write the EGF of 2ⁿ. Use EGFs to show that the number of strings of length n over {a, b, c} that contain at least one a is 3ⁿ − 2ⁿ.
2. (Intermediate) Show that the number of strings of length n over {0, 1, 2} with an even number of 0s and an odd number of 1s is (3ⁿ − (−1)ⁿ)/4.
3. (Interview) In how many ways can n people be split into teams of at least 2 people each, where the teams are unlabeled? Give the EGF and a recurrence you could turn into an O(n²) dynamic program. Hint: a team is a set of size at least 2, with EGF eˣ − 1 − x, so the answer has EGF e^{eˣ − 1 − x}. For the recurrence, follow the Bell-number recurrence of 5.6, but require the team of person n + 1 to contain at least one other person.

<a id="7-discrete-probability"></a>