Integer Partitions, in depth: Ferrers diagrams, the O(n^2) DP, Euler's pentagonal recurrence and generating partitions

Ferrers diagrams and conjugation

Write a partition with parts in non-increasing order and draw part i as a row of i dots: the Ferrers diagram. Reading the columns instead of the rows gives the conjugate partition. Because conjugation is its own inverse, it is a bijection, and it swaps two properties: the number of parts becomes the largest part. That single picture proves that partitions of n into at most k parts are equinumerous with partitions of n into parts no larger than k, which matters below because the second form is the one with the easy recurrence.

Counting with a part-size DP

Let P(m, k) be the number of partitions of m using parts no larger than k. Either the partition uses no part equal to k, giving P(m, k - 1), or it uses at least one, and removing one copy leaves a partition of m - k with parts no larger than k, giving P(m - k, k). Rolling the k dimension away gives a loop identical to counting coin-change combinations with coins 1, 2, ..., n:

def partition_counts(n):
    # p[m] = number of partitions of m, for every m in 0..n. O(n^2) time, O(n) space.
    p = [1] + [0] * n                 # one way to partition 0: the empty partition
    for k in range(1, n + 1):         # allow part size k
        for m in range(k, n + 1):
            p[m] += p[m - k]
    return p

assert partition_counts(10)[10] == 42

The outer loop must be over part sizes and the inner over totals. Swap them and each total absorbs every ordering of its parts, so you count compositions instead, the same bug that turns combinations into permutations in coin change. Restricting the parts is just a matter of which k the outer loop visits: only odd k, only parts from a set, or k up to some bound. Keeping the full table instead gives P(m, k) for every k; for example the partitions of 10 into at most 3 parts number P(10, 3) = 14. That table is also what you need for exact uniform sampling, below.

Trace it for n = 5 to see the mechanics. The array starts as [1, 0, 0, 0, 0, 0]. Allowing parts of size 1 fills it with ones: there is exactly one partition of each m into 1s. Allowing size 2 adds p[m - 2] to each entry from m = 2 upwards, giving [1, 1, 2, 2, 3, 3]; for instance 4 now counts 1+1+1+1, 2+1+1 and 2+2. Size 3 gives [1, 1, 2, 3, 4, 5], size 4 gives [1, 1, 2, 3, 5, 6], and size 5 adds the single partition 5 itself, ending at [1, 1, 2, 3, 5, 7]. Each row after pass k is the column P(m, k) of the full table, which is why keeping every pass costs O(n^2) memory but gives the restricted counts for free.

Euler's pentagonal recurrence

When n is large the O(n^2) loop is wasteful. Euler's pentagonal number theorem says that the product of (1 - x^k) over all k has almost every coefficient zero: it equals the sum over all integers j of (-1)^j x^(j(3j - 1)/2). Since the generating function of p(n) is the reciprocal of that product, multiplying the two gives a recurrence with only O(sqrt n) terms:

p(m) = sum_{j >= 1} (-1)^(j+1) * [ p(m - j(3j-1)/2) + p(m - j(3j+1)/2) ]
#      generalised pentagonal offsets: 1, 2, 5, 7, 12, 15, 22, 26, ...
#      signs, in pairs:                +  +  -  -  +   +   -   -

def partition_counts_mod(n, mod=10**9 + 7):
    offsets = []                                  # (offset, sign) in increasing order
    j = 1
    while j * (3 * j - 1) // 2 <= n:
        sign = 1 if j % 2 else -1
        offsets.append((j * (3 * j - 1) // 2, sign))
        offsets.append((j * (3 * j + 1) // 2, sign))
        j += 1
    p = [1] + [0] * n
    for m in range(1, n + 1):
        total = 0
        for off, sign in offsets:
            if off > m:
                break
            total += sign * p[m - off]
        p[m] = total % mod
    return p

There are 516 generalised pentagonal numbers up to 100,000, so each p(m) costs at most 516 lookups and the full table is about 3.4 x 10^7 operations: fine in C, slow in pure Python (one run took about half a minute). For p(n) modulo a prime and every n up to N at once, the fastest route is to invert the power series of the product of (1 - x^k), which is sparse thanks to the theorem, using NTT-based multiplication in O(N log N). The sign pattern is the usual source of bugs: two plus, two minus, repeating. Always cross-check against the O(n^2) loop for small n; the two agree for every n up to 1,000.

How fast p(n) grows

Hardy and Ramanujan showed in 1918 that p(n) is asymptotic to exp(pi sqrt(2n/3)) / (4 n sqrt 3). The ratio of the approximation to the true value is 1.145 at n = 10, 1.046 at n = 100 and 1.014 at n = 1000, so it is a good size estimate but never an exact answer; Rademacher's convergent series gives the exact value, and is what libraries use to compute a single p(n) for huge n. Practically: p(n) first exceeds a signed 64-bit integer at n = 406 and an unsigned one at n = 417. Use big integers, a modulus or floating-point logs beyond that, and add an assertion so a silent wrap cannot hide.

Ramanujan also found that p(5k + 4) is divisible by 5, p(7k + 5) by 7 and p(11k + 6) by 11 for every k. They make excellent unit tests: running them over the first thousand values catches almost any off-by-one in a recurrence.

Generating and sampling partitions

Sometimes you need the partitions themselves, for example to enumerate every way to split a batch of identical shards among unlabelled workers. Generate them in reverse lexicographic order: start from [n]; at each step strip the trailing 1s, decrease the last remaining part by one, and refill with copies of that new value plus a remainder.

def partitions(n):
    a = [n]
    while True:
        yield tuple(a)
        rem = 0
        while a and a[-1] == 1:       # collect trailing ones
            a.pop()
            rem += 1
        if not a:
            return                    # we just yielded 1+1+...+1
        a[-1] -= 1
        rem += 1
        k = a[-1]
        while rem > k:                # refill with parts of size k
            a.append(k)
            rem -= k
        if rem:
            a.append(rem)

list(partitions(5))
# [(5,), (4, 1), (3, 2), (3, 1, 1), (2, 2, 1), (2, 1, 1, 1), (1, 1, 1, 1, 1)]

Because it pops the trailing 1s one at a time, this simple version does about sqrt(n) work per partition on average; Knuth's Algorithm P and the ZS1 algorithm reach constant amortised time by tracking the last part larger than 1 instead of storing the 1s. Its count matches the formula: it yields 5,604 partitions for n = 30, which is p(30). Remember the growth table before enumerating: p(100) is 190 million, so anything beyond n of about 80 should count, sample or prune rather than list. To sample a partition uniformly at random, use the P(m, k) table: choose the largest part k with probability proportional to the number of completions, subtract it, and repeat with parts no larger than k. That is the same unranking idea as for combinations in the combinatorics article.

Restricted partitions and generating functions

Generating functions organise all the variants. Partitions of n are the coefficients of the product over k of 1 / (1 - x^k); each factor chooses how many copies of k to use. Restrict the factors and you restrict the parts. Distinct parts give the product of (1 + x^k); odd parts give the product of 1 / (1 - x^k) over odd k. Euler noticed the two products are equal, because (1 + x^k) = (1 - x^2k) / (1 - x^k) and the even factors cancel, so the number of partitions into distinct parts always equals the number into odd parts. For n = 1 to 10 both counts run 1, 1, 2, 2, 3, 4, 5, 6, 8, 10. In code, distinct parts is the 0/1 version of the DP loop (iterate totals downwards so each part is used at most once), and odd parts is the unbounded loop over odd k only. Checking one against the other is another cheap test.

Failure modes

• Counting compositions by accident. Loop order in the DP decides whether 1+2 and 2+1 are the same. Test p(4) = 5, not 8.
• Overflow. 64-bit integers fail at n = 406. Languages without big integers wrap silently.
• Pentagonal sign errors. A wrong sign pattern gives plausible-looking numbers that drift from n = 5 onwards. Compare with the DP for every n up to a few hundred.
• Enumerating when you should count. Listing partitions of 120 means about 1.8 billion tuples. Estimate with the growth table first.
• Recursion depth. A naive recursive P(m, k) without memoisation is exponential, and with memoisation in Python it hits the recursion limit near n = 1000. Use the iterative table.
• Confusing the problems. "Partition" also names graph partitioning, set partitions (Bell numbers) and the NP-complete equal-sum split. Check which one a spec means.

Trade-offs

What to do next

1. Implement partition_counts and assert p(4) = 5, p(10) = 42 and p(100) = 190,569,292.
2. Add the pentagonal version and check that it agrees with the DP up to n = 1,000, plus the three Ramanujan congruences.
3. Write the distinct-parts and odd-parts counters and check they agree for n up to 200.
4. Use the generator on n = 8 by hand, then sample uniformly from the P(m, k) table and histogram the results against the exact counts.
5. Decide the number type before shipping: big integers, a modulus, or a hard cap at n = 405.
6. Read the combinatorics article for the other counting families, then the NTT article if you need p(n) modulo a prime for very large n.