Catalan Numbers, in depth: the first-return recurrence, the reflection proof, bijections, exact and modular computation, ranking and uniform random trees

The first-return recurrence

Take a balanced string of n bracket pairs, or equivalently a Dyck path of n up steps and n down steps that never dips below the axis. Its first character is an opening bracket. Find the bracket that closes it. Between them sits a balanced string of some size i, and after the closer sits a balanced string of size n - 1 - i. The split is unique and every pair of smaller strings produces exactly one larger string, so

C₀ = 1, and C_n = sum over i from 0 to n-1 of C_i C_n-1-i.

Check it by hand: C₃ = C₀C₂ + C₁C₁ + C₂C₀ = 2 + 1 + 2 = 5, and the five strings are ((())), (()()), (())(), ()(()) and ()()(). This recurrence is the definition that transfers to other structures: any object that splits uniquely into an ordered pair of smaller objects of the same kind, with total size reduced by one, is counted by the Catalan numbers.

The closed form by reflection

The closed form is C_n = C(2n, n) / (n + 1). The cleanest proof counts bad paths. A path of n ups and n downs, ignoring the floor, is a choice of which n of 2n steps go up: C(2n, n) paths. A bad path touches height -1 somewhere. Take the first such touch and flip every step after it, ups to downs and downs to ups. The path that ended at height 0 now ends at height -2, so it has n - 1 ups and n + 1 downs. The flip is reversible, because any path ending at -2 must touch -1, and flipping after that first touch restores the original. So bad paths number C(2n, n + 1), and

C_n = C(2n, n) - C(2n, n+1) = C(2n, n) / (n + 1).

The same argument generalises. Paths with a ups and b downs, a at least b, that never go below zero number C(a+b, b) - C(a+b, b-1). That is the ballot problem, and it is the tool for counting balanced prefixes, which the ranking code below needs.

One count, many structures

Each structure below satisfies the first-return recurrence, which is why the counts agree. Knowing the bijection lets you convert a problem you cannot count into one you can.

The binary-tree row is the one engineers meet most: the number of distinct binary search tree shapes for n keys is C_n, which is why a search over all BSTs, as opposed to the dynamic program in optimal binary search trees, is hopeless. The parenthesisation row is why a product of n matrices has C_n-1 evaluation orders, the space that matrix chain multiplication searches in polynomial time instead of enumerating.

Computing the numbers

There are three ways to compute C_n, and the right one depends on how big n is and what arithmetic you have.

def catalan_dp(n):                       # O(n^2), exactly the recurrence
    C = [0] * (n + 1)
    C[0] = 1
    for m in range(1, n + 1):
        C[m] = sum(C[i] * C[m - 1 - i] for i in range(m))
    return C

def catalan_exact(n):                    # O(n) big-integer multiplications
    c = 1
    for k in range(n):                   # C(k+1) = C(k) * 2(2k+1) / (k+2)
        c = c * 2 * (2 * k + 1) // (k + 2)
    return c

def catalan_mod(n, p=10**9 + 7):         # O(n), for prime p > 2n
    fact = [1] * (2 * n + 2)
    for i in range(1, 2 * n + 2):
        fact[i] = fact[i - 1] * i % p
    inv = lambda x: pow(x, p - 2, p)     # Fermat inverse
    return fact[2 * n] * inv(fact[n]) % p * inv(fact[n + 1]) % p

The multiplicative step is exact because the division always comes out whole: C_k+1 is an integer, and the code multiplies first and divides second. Reverse that order and you truncate. Running these: catalan_dp(15) ends 742900, 2674440, 9694845; catalan_exact(35) is 3,116,285,494,907,301,262; and catalan_mod(1000) returns 110961515, which matches the exact value reduced modulo 10⁹ + 7.

64-bit limits. C₃₅ is the largest value that fits a signed 64-bit integer. C₃₆ = 11,959,798,385,860,453,492 fits only unsigned, and C₃₇ fits neither. But the multiplicative loop fails earlier: computing C₃₄ from C₃₃ = 212,336,130,412,243,110 requires the intermediate product C₃₃ times 134, about 2.8 times 10¹⁹, which overflows signed 64-bit even though C₃₄ itself fits. In a fixed-width language, either use 128-bit intermediates or divide by the greatest common divisor before multiplying; or simply precompute the 36 values into a table.

For catalan_mod, the prime must exceed 2n so that no factorial is zero modulo p. For smaller primes you need Lucas-style digit decomposition, which our combinatorics article covers.

Generating, ranking and unranking

To enumerate every object, generate balanced strings by backtracking with two counters: you may open while fewer than n are open, and close while closed is less than opened. Every leaf of the search is valid, so the cost is proportional to the output, C_n strings of length 2n. The same pruning pattern appears in backtracking over permutations and combinations.

def balanced(n):
    out = []
    def go(prefix, opened, closed):
        if len(prefix) == 2 * n:
            out.append("".join(prefix)); return
        if opened < n:
            prefix.append("("); go(prefix, opened + 1, closed); prefix.pop()
        if closed < opened:
            prefix.append(")"); go(prefix, opened, closed + 1); prefix.pop()
    go([], 0, 0)
    return out

# balanced(3) -> ['((()))', '(()())', '(())()', '()(())', '()()()']

Ranking and unranking. Often you do not want all of them, only the k-th, for sharding a search across machines, or for using an integer as a compact key. Precompute W[r][h], the number of ways to finish with r symbols left from height h without dipping below zero. Then walk the string: wherever it has a closer, every word that has an opener at that position comes earlier in lexicographic order, so add their count.

def completions(n):
    W = [[0] * (2 * n + 2) for _ in range(2 * n + 1)]
    W[0][0] = 1
    for r in range(1, 2 * n + 1):
        for h in range(0, 2 * n + 1):
            W[r][h] = W[r - 1][h + 1] + (W[r - 1][h - 1] if h > 0 else 0)
    return W

def rank(word, W):
    r, h = 0, 0
    for i, ch in enumerate(word):
        left = len(word) - i - 1
        if ch == ")":
            r += W[left][h + 1]          # words with "(" here sort first
            h -= 1
        else:
            h += 1
    return r

def unrank(r, n, W):
    out, h = [], 0
    for i in range(2 * n):
        cnt = W[2 * n - i - 1][h + 1]
        if r < cnt:
            out.append("("); h += 1
        else:
            r -= cnt; out.append(")"); h -= 1
    return "".join(out)

Tested at n = 4, rank maps the 14 strings from balanced(4) to 0 through 13 in order and unrank inverts it exactly; W[20][0] for n = 10 is 16796, which is C₁₀, a useful self-check. Both run in O(n) after an O(n²) table.

Uniform random objects

Random testing of a parser, a tree library or a query planner needs uniformly random trees, and the obvious approach, flipping a coin at each step and rejecting failures, is biased toward shallow shapes. The cycle lemma gives an exact method. Take n up steps and n + 1 down steps in uniformly random order; the total is -1. Of the 2n + 1 rotations of that sequence, exactly one has every proper prefix sum non-negative: the one that starts just after the first position where the running sum hits its minimum. Drop that rotation's final down step and you have a uniformly random Dyck word.

import random

def random_dyck(n, rng=random):
    steps = [1] * n + [-1] * (n + 1)
    rng.shuffle(steps)
    s, low, cut = 0, 0, 0
    for i, x in enumerate(steps):
        s += x
        if s < low:                      # first time the minimum is reached
            low, cut = s, i + 1
    rotated = steps[cut:] + steps[:cut]
    return "".join("(" if x == 1 else ")" for x in rotated[:-1])

Uniformity follows because each of the C(2n+1, n) sequences maps to one Dyck word and every Dyck word receives exactly 2n + 1 of them. Empirically, 140,000 draws at n = 4 with a fixed seed hit all 14 words, each between 9,864 and 10,187 times against an expected 10,000. Convert the word to a binary tree with the first-return split and you have a uniform random tree in O(n).

Growth and why enumeration dies

Asymptotically C_n is about 4ⁿ / (n^3/2 sqrt(pi)). The estimate overshoots, by 11 percent at n = 10, 6 percent at n = 20 and 2 percent at n = 50, but the shape is what matters: each extra element multiplies the count by nearly four.

That table is the practical lesson. Any algorithm that enumerates bracketings, tree shapes or evaluation orders is exponential with base four. When the objective decomposes along the first-return split, which it does whenever the cost of a tree is the cost of its subtrees plus a root term, the same recurrence becomes an interval dynamic program over O(n²) subproblems; see dynamic programming for the general pattern.

Failure modes

• Divide before multiply. Writing c = c / (k + 2) * 2 * (2k + 1) truncates: the very first step computes 1 // 2 = 0, and every value after that is 0 in integer arithmetic.
• Silent 64-bit overflow. The loop goes wrong at C₃₄, not C₃₆; tests that stop at n = 20 never see it.
• Modulus too small. catalan_mod with p not greater than 2n hits a zero factorial and returns 0 without complaint.
• Off-by-one family. n matrices give C_n-1 orders and an n-gon has C_n-2 triangulations; index from the recurrence, not from memory.
• Biased random trees. Growing trees by random insertion or coin flips is not uniform over shapes; use the cycle lemma.
• Using the asymptotic as a value. It is an estimate with percent-level error at moderate n; never use it where an exact count is needed.

Trade-offs

What to do next

1. Implement catalan_dp and catalan_exact and assert they agree for n up to 30.
2. Port the multiplicative loop to a 64-bit language and confirm it breaks at C₃₄; fix it with 128-bit intermediates.
3. Write balanced(n) and check the lengths of its output against your table for n up to 10.
4. Build rank and unrank and use them to shard an enumeration across workers by integer range.
5. Use random_dyck to generate uniform random binary trees and fuzz a tree or parser routine with them.
6. When you next see a search over bracketings or tree shapes, write the first-return split and turn it into an interval DP.