14 — Dynamic programming

Dynamic programming (DP) solves a problem by splitting it into smaller versions of itself, solving each smaller version once, and storing the answer so it never has to be computed again. It’s the technique behind spell checkers (edit distance), diff (longest common subsequence), sequence alignment in bioinformatics, the Viterbi decoder, the PostgreSQL planner’s search for a join order, and exact solvers for small traveling salesman instances.

The name was coined by Richard Bellman in the 1950s. “Programming” meant planning, as in “linear programming”, not writing code. “Dynamic” referred to decisions made in stages over time. Neither word helps much with understanding the technique, which comes down to two properties of a problem and a disciplined way of exploiting them.

Timings in this chapter were measured with Python 3.14 on one Mac. Exact numbers will differ elsewhere, but the ratios between approaches are what matter. Every implementation was checked against a brute-force solution on hundreds of random small inputs.

What makes a problem a DP problem

A problem is a good fit for DP when it has both of these properties.

Optimal substructure. An optimal solution to the whole problem is built from optimal solutions to subproblems. The shortest path from A to C through B consists of a shortest path from A to B followed by a shortest path from B to C. If it didn’t, swapping in a shorter A-to-B segment would give a shorter A-to-C path, a contradiction. That “cut and paste” argument is how optimal substructure is usually proved.

Not every problem has it. The longest simple path (no repeated vertices) from A to C through B is not the longest simple path from A to B plus the longest from B to C, because the two halves can share vertices and the combination is no longer simple. The subproblems interfere with each other. That’s why longest simple path is NP-hard while shortest path is polynomial.

Overlapping subproblems. A naive recursive solution asks the same subquestion many times. Fibonacci is the standard example:

                          fib(5)
                 /                    \
            fib(4)                    fib(3)
           /      \                  /      \
       fib(3)    fib(2)          fib(2)    fib(1)
       /    \    /    \          /    \
   fib(2) fib(1) f(1) f(0)    f(1)   f(0)
   /    \
 f(1)  f(0)

fib(3) is computed twice, fib(2) three times. The number of calls grows like $\varphi^n$ where $\varphi \approx 1.618$, yet there are only $n + 1$ distinct subproblems. Storing each answer the first time collapses the exponential tree into a linear chain.

If subproblems don’t overlap, storing answers gains nothing. Merge sort has optimal substructure but every subarray is sorted exactly once, so it’s plain divide and conquer, not DP.

DP is also different from greedy algorithms (chapter 16). A greedy algorithm commits to one locally best choice at each step and never revisits it. DP evaluates every choice at each step and keeps the best, relying on the stored subproblem answers to make that affordable. When a greedy choice is provably safe, greedy is faster. When it isn’t, DP is the fallback that still runs in polynomial time.

The recipe

Almost every DP solution is designed by answering five questions in order.

1. State. What does a subproblem look like, and which parameters identify it? Write it as a sentence: “dp[i][c] is the best value achievable using only the first i items with capacity c.” If that sentence can’t be written precisely, nothing after it will work.
2. Transition. How is a state’s answer computed from smaller states? This usually comes from considering the last decision: was the last item taken or not, did the last characters match or not, where was the last cut made.
3. Base cases. Which states have answers that don’t depend on anything smaller? Usually the empty prefix, a zero capacity, or a single element.
4. Order. In what order must the states be computed so that every state’s dependencies are ready when it’s evaluated? Recursion with memoization figures this out automatically. A table needs it spelled out.
5. Answer. Which state, or which combination of states, holds the final result? Often dp[n] or dp[n][m], but sometimes it’s the maximum over all states, as in longest increasing subsequence.

The cost of a DP is almost always

$\text{time} = (\text{number of states}) \times (\text{work per transition})$

and the space is the number of states kept alive at once. That formula is the fastest way to tell whether a state design is viable before writing any code: a state with three parameters each up to 1,000 means a billion states, which no amount of tuning will make fast in Python.

A 1D example: House Robber

A row of houses holds amounts of money. Robbing two adjacent houses triggers an alarm. Find the largest total that can be taken.

• State: dp[i] is the best total from the first i houses.
• Transition: house i - 1 is either skipped, giving dp[i - 1], or robbed, in which case house i - 2 must be skipped, giving dp[i - 2] + houses[i - 1]. Take the maximum.
• Base cases: dp[0] = 0 (no houses), and dp[1] = houses[0] follows from the transition when dp[-1] is treated as 0.
• Order: increasing i.
• Answer: dp[n].

def rob_with_choice(houses):
    n = len(houses)
    dp = [0] * (n + 1)          # dp[i] = best using the first i houses
    for i in range(1, n + 1):
        dp[i] = max(dp[i - 1], (dp[i - 2] if i >= 2 else 0) + houses[i - 1])
    picked, i = [], n
    while i > 0:
        if dp[i] == dp[i - 1]:
            i -= 1                     # house i-1 was not taken
        else:
            picked.append(i - 1)       # house i-1 was taken, skip its neighbor
            i -= 2
    return dp[n], picked[::-1]

>>> rob_with_choice([2, 7, 9, 3, 1])
(12, [0, 2, 4])

The second half of the function is reconstruction, covered below. The value alone needs only the last two entries of dp, which leads to the constant-space version:

def rob(houses):
    take, skip = 0, 0          # best ending with house i taken / not taken
    for x in houses:
        take, skip = skip + x, max(take, skip)
    return max(take, skip)

Time is $O(n)$ for both. Space is $O(n)$ for the table version and $O(1)$ for the rolling version.

Memoization and tabulation

There are two ways to turn a recurrence into an algorithm.

Memoization (top-down) keeps the natural recursive function and adds a cache. The first call for each argument computes and stores the result. Later calls return the stored value. In Python, functools.cache (an unbounded lru_cache) does this in one line:

from functools import cache

def fib_naive(n):
    if n < 2:
        return n
    return fib_naive(n - 1) + fib_naive(n - 2)

@cache
def fib_memo(n):
    if n < 2:
        return n
    return fib_memo(n - 1) + fib_memo(n - 2)

Tabulation (bottom-up) allocates a table and fills it in an order that guarantees each entry’s dependencies are already filled:

def fib_table(n):
    if n < 2:
        return n
    dp = [0] * (n + 1)
    dp[1] = 1
    for i in range(2, n + 1):
        dp[i] = dp[i - 1] + dp[i - 2]
    return dp[n]

Both compute each subproblem once, so they have the same asymptotic time. They differ in practice:

Memoization wins when only a sparse subset of the state space is reachable, for example a state like (position, remaining_budget) where most budget values never occur. Tabulation wins on raw speed when most states are needed anyway, and it’s the only one of the two that can drop old rows to save memory.

Measured: naive, memoized and tabulated

Fibonacci, with the cache cleared before each memoized run:

The naive column was not run at 500 and 900: extrapolating the growth rate, fib_naive(500) would make more calls than there are atoms in the universe.

The naive version multiplies by about $\varphi^5 \approx 11$ every time n grows by 5, exactly as the call tree predicts. Memoization and tabulation are both linear, with the table several times faster because a list index is much cheaper than a Python function call plus a dictionary lookup.

The same comparison on a 2D problem, the length of the longest common subsequence of two random DNA strings (the recurrence is explained in the LCS section):

For the short rows, the two strings were drawn from disjoint alphabets (A/C against G/T), so no characters ever match and the naive recursion branches at every step, its worst case. On random DNA strings with matches the naive version was faster (3.28 ms at length 12) but still exponential.

Naive recursion is exponential: two more characters per string and the time grows by a factor of about 14. Both DP versions are $O(mn)$. At 1,000 by 1,000 there are a million states and the memoized version needed the recursion limit raised to even run.

The recursion limit

CPython limits recursion depth to protect the C stack. The default is 1,000:

>>> import sys
>>> sys.getrecursionlimit()
1000
>>> fib_memo(5000)
Traceback (most recent call last):
  ...
RecursionError: maximum recursion depth exceeded

A memoized recursion that goes n levels deep before hitting a base case needs n stack frames, whatever the cache holds. Three ways around it:

• Call sys.setrecursionlimit(...) with a larger value. This works for moderate depths, but very deep recursion can still crash the interpreter with a real stack overflow, which is worse than an exception.
• Warm the cache from the bottom: call fib_memo(500), fib_memo(1000), and so on in increasing order, so no single call recurses far. Warming the cache in steps of 500 was enough to reach fib_memo(5000), a 1,045-digit number, under the default limit.
• Rewrite it as a table. This is the robust answer for any DP whose depth scales with input size.

Space optimization with rolling arrays

When row i of a table depends only on row i - 1 (or on a constant number of previous rows), the older rows can be discarded. Fibonacci needs two numbers:

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

For a 2D table, two rows suffice, and sometimes one row updated in place. LCS length with two rows, swapping the strings so that the row is as short as possible:

def lcs_length_two_rows(a, b):
    if len(b) > len(a):
        a, b = b, a                       # keep the row as short as possible
    prev = [0] * (len(b) + 1)
    for ch in a:
        cur = [0] * (len(b) + 1)
        for j in range(1, len(b) + 1):
            if ch == b[j - 1]:
                cur[j] = prev[j - 1] + 1
            else:
                cur[j] = max(prev[j], cur[j - 1])
        prev = cur
    return prev[-1]

Space falls from $O(mn)$ to $O(\min(m, n))$. The catch is that the discarded rows are exactly what reconstruction walks back through. If the actual solution is needed, either keep the full table, store a compact “which choice was made” table (one byte per state instead of a number), or use Hirschberg’s divide-and-conquer method, which recovers an LCS or alignment in $O(mn)$ time and linear space by computing the middle row from both ends and recursing on the two halves.

Reconstructing the solution, not just its value

The table records the best value for every state. To recover the choices that produced the final answer, start at the answer state and ask, at each step, which transition produced the stored value. Move to that predecessor, record the choice, and repeat until a base case is reached. The choices come out in reverse order.

Two ways to implement it:

• Recompute the decision by comparing the stored value with each candidate, as rob_with_choice does above (dp[i] == dp[i - 1] means house i - 1 was skipped). This needs no extra memory.
• Store a parent pointer or choice for each state during the forward pass, as in the LIS and Held-Karp code below. This is more direct when the transition has many candidates or when ties make recomputation ambiguous.

When several optimal solutions exist, the tie-breaking rule in the walk decides which one is returned. That’s fine as long as the walk only ever follows transitions that achieve the stored value.

1D, 2D and 3D DP

The dimension of a DP is the number of parameters in its state.

Higher dimensions are not harder in principle, only more expensive. Each added parameter multiplies the number of states. A 3D example, the longest common subsequence of three strings, extends the two-string recurrence by one index:

def lcs3(a, b, c):
    dp = [[[0] * (len(c) + 1) for _ in range(len(b) + 1)] for _ in range(len(a) + 1)]
    for i in range(1, len(a) + 1):
        for j in range(1, len(b) + 1):
            for k in range(1, len(c) + 1):
                if a[i - 1] == b[j - 1] == c[k - 1]:
                    dp[i][j][k] = dp[i - 1][j - 1][k - 1] + 1
                else:
                    dp[i][j][k] = max(dp[i - 1][j][k], dp[i][j - 1][k], dp[i][j][k - 1])
    return dp[-1][-1][-1]

>>> lcs3("AGGT12", "12TXAYB", "12XBA")
2

It runs in $O(\lvert a \rvert \lvert b \rvert \lvert c \rvert)$ time and space. One subtlety: the LCS of three strings is not the LCS of the first two intersected with the third. Computing a two-string LCS first and then matching it against the third string can miss the optimum, because the first step may pick the wrong one of several equally long common subsequences. With BBA, ABAB and AABAA, the two-string lcs above returns BB for the first pair, which shares only one character with the third string, while BA is common to all three.

The knapsack family

A large group of DP problems share one shape: a list of items, a numeric budget, and for each item a decision about whether (or how many times) to use it. The budget becomes one dimension of the state, which makes these algorithms pseudo-polynomial: polynomial in the numeric value of the budget, not in the number of bits needed to write it down. A capacity of $10^{12}$ is a 13-digit input but a trillion-column table. The underlying decision problems (knapsack, subset sum, partition) are NP-complete, and these DPs don’t contradict that, because their running time is exponential in the input’s bit length.

0/1 knapsack

Each of n items has a weight and a value. Choose a subset with total weight at most capacity and maximum total value. Each item can be taken at most once.

• State: dp[i][c] is the best value using only the first i items with capacity c.
• Transition: item i - 1 is either skipped, giving dp[i - 1][c], or taken (if it fits), giving dp[i - 1][c - w] + v.
• Base case: dp[0][c] = 0 for every c.
• Answer: dp[n][capacity].

def knapsack_2d(weights, values, capacity):
    n = len(weights)
    dp = [[0] * (capacity + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        w, v = weights[i - 1], values[i - 1]
        for c in range(capacity + 1):
            dp[i][c] = dp[i - 1][c]                       # skip item i-1
            if w <= c:
                dp[i][c] = max(dp[i][c], dp[i - 1][c - w] + v)   # take it
    chosen, c = [], capacity
    for i in range(n, 0, -1):
        if dp[i][c] != dp[i - 1][c]:       # value changed, so item i-1 was taken
            chosen.append(i - 1)
            c -= weights[i - 1]
    return dp[n][capacity], chosen[::-1]

With weights [1, 3, 4, 5], values [1, 4, 5, 7] and capacity 7, the table is:

capacity →   0   1   2   3   4   5   6   7
no items     0   0   0   0   0   0   0   0
+ item 0     0   1   1   1   1   1   1   1     (w=1, v=1)
+ item 1     0   1   1   4   5   5   5   5     (w=3, v=4)
+ item 2     0   1   1   4   5   6   6   9     (w=4, v=5)
+ item 3     0   1   1   4   5   7   8   9     (w=5, v=7)

>>> knapsack_2d([1, 3, 4, 5], [1, 4, 5, 7], 7)
(9, [1, 2])

The walk back starts at the bottom-right 9. The row above also holds 9, so item 3 wasn’t taken. Moving up, 9 differs from the 5 above it, so item 2 was taken and capacity drops to 3. At row ”+ item 1” column 3, the value 4 differs from the 1 above, so item 1 was taken and capacity drops to 0. Items 1 and 2 weigh 7 and are worth 9. A greedy choice by value-to-weight ratio would take item 3 (ratio 1.4) first, then item 0, for a total of 8, which is why knapsack needs DP.

Time is $O(n \cdot C)$ for capacity $C$, and space is $O(n \cdot C)$ with reconstruction.

Why the 1D version loops backwards

Row i reads only row i - 1, so one array can hold both: before an entry is overwritten it holds the previous row’s value. The order of the inner loop then decides which row each read sees.

def knapsack_1d(weights, values, capacity):
    dp = [0] * (capacity + 1)
    for w, v in zip(weights, values):
        for c in range(capacity, w - 1, -1):      # backwards: each item used once
            dp[c] = max(dp[c], dp[c - w] + v)
    return dp[capacity]

The update at c reads dp[c - w], an entry to its left. Going from high c to low, that entry hasn’t been updated yet in this pass, so it still holds the “without this item” value from the previous row. Going from low to high, dp[c - w] has already been updated and may already include the current item, so the item can be added again.

A single item with weight 2 and value 3, capacity 6:

backwards  [0, 0, 3, 3, 3, 3, 3]    item used at most once: best is 3
forwards   [0, 0, 3, 3, 6, 6, 9]    dp[4] reused dp[2], which already had the item

def knapsack_1d_forwards_bug(weights, values, capacity):
    dp = [0] * (capacity + 1)
    for w, v in zip(weights, values):
        for c in range(w, capacity + 1):          # forwards: item can be reused
            dp[c] = max(dp[c], dp[c - w] + v)
    return dp[capacity]

>>> knapsack_1d([2], [3], 6), knapsack_1d_forwards_bug([2], [3], 6)
(3, 9)

The forward loop doesn’t produce garbage, it produces the answer to a different problem: the unbounded knapsack, where each item can be used any number of times. The verification run confirmed that the forward 1D loop matches an independent unbounded-knapsack solver on every random input. The 1D version uses $O(C)$ space but loses the ability to reconstruct the chosen items, since earlier rows are gone.

Unbounded knapsack

With unlimited copies of each item, the state no longer needs the item index. dp[c] is the best value with capacity c, and the last item added can be any item that fits:

def unbounded_knapsack(weights, values, capacity):
    dp = [0] * (capacity + 1)
    for c in range(1, capacity + 1):
        for w, v in zip(weights, values):
            if w <= c:
                dp[c] = max(dp[c], dp[c - w] + v)
    return dp[capacity]

Time $O(n \cdot C)$, space $O(C)$. The forward 1D loop from the previous section computes the same values.

Subset sum and partition equal subset sum

Subset sum asks whether some subset of the numbers adds up to exactly target. It’s a 0/1 knapsack where the value is a boolean “reachable”, so the inner loop again runs backwards:

def subset_sum(nums, target):
    reachable = [False] * (target + 1)
    reachable[0] = True
    for x in nums:
        for s in range(target, x - 1, -1):
            if reachable[s - x]:
                reachable[s] = True
    return reachable[target]

In Python, an arbitrary-precision integer makes a fast bitset. Bit s set means sum s is reachable, and adding a number x to every reachable sum is a shift. Shifting a whole integer processes all sums at once in C. On 200 random numbers up to 1,000 with a target of 50,061, the list version took 165 ms and the bitset version 326 µs:

def subset_sum_bits(nums, target):
    mask = 1                        # bit s set means sum s is reachable
    for x in nums:
        mask |= mask << x
    return (mask >> target) & 1 == 1

Partition Equal Subset Sum asks whether the numbers can be split into two groups with equal sums. That’s possible only if the total is even and some subset sums to half of it:

def can_partition(nums):
    total = sum(nums)
    if total % 2:
        return False
    return subset_sum(nums, total // 2)

>>> can_partition([1, 5, 11, 5]), can_partition([1, 2, 3, 5])
(True, False)

The first splits as [1, 5, 5] and [11]. The second has an odd total of 11.

Target sum

Target Sum puts a + or - in front of each number and counts the sign assignments that evaluate to target. It looks like a search over $2^n$ sign patterns, but a little algebra turns it into subset counting. Let $P$ be the sum of the numbers given + and $N$ the sum of those given -. Then $P - N = \text{target}$ and $P + N = \text{total}$, so $P = (\text{total} + \text{target}) / 2$. The question becomes: how many subsets sum to $P$?

def target_sum_ways(nums, target):
    total = sum(nums)
    if abs(target) > total or (total + target) % 2:
        return 0
    goal = (total + target) // 2          # sum of the numbers given a + sign
    ways = [0] * (goal + 1)
    ways[0] = 1
    for x in nums:
        for s in range(goal, x - 1, -1):
            ways[s] += ways[s - x]
    return ways[goal]

>>> target_sum_ways([1, 1, 1, 1, 1], 3)
5

Zeros need care. A zero can take either sign without changing the sum, so it doubles the count. The loop above handles this correctly (with x = 0 the inner loop adds ways[s] to itself), and the random tests included zeros. A common bug is a loop bound like range(goal, x, -1) that skips s = x and silently drops those doublings.

Coin change: fewest coins

Given coin denominations with unlimited supply, find the fewest coins that make amount, or report that it can’t be done. It’s an unbounded knapsack minimizing count:

from math import inf

def coin_change_min(coins, amount):
    dp = [0] + [inf] * amount
    for a in range(1, amount + 1):
        for coin in coins:
            if coin <= a and dp[a - coin] + 1 < dp[a]:
                dp[a] = dp[a - coin] + 1
    return dp[amount] if dp[amount] != inf else -1

>>> coin_change_min([1, 3, 4], 6)
2
>>> coin_change_min([2], 3)
-1

With coins [1, 3, 4], greedy takes 4, then 1, then 1, using three coins. The optimum is 3 + 3. The greedy method is optimal for some coin systems (called canonical, which includes most real currencies) but not in general. Time $O(n \cdot \text{amount})$, space $O(\text{amount})$. Initializing unreachable amounts to infinity, not to 0 or -1, keeps the min correct without special cases.

Coin change: number of ways, and why loop order matters

Counting the ways to make amount has two versions that differ only in the order of the two loops.

def coin_combinations(coins, amount):
    ways = [1] + [0] * amount
    for coin in coins:                      # coins outer: each multiset once
        for a in range(coin, amount + 1):
            ways[a] += ways[a - coin]
    return ways[amount]


def coin_permutations(coins, amount):
    ways = [1] + [0] * amount
    for a in range(1, amount + 1):          # amounts outer: every order counted
        for coin in coins:
            if coin <= a:
                ways[a] += ways[a - coin]
    return ways[amount]

>>> coin_combinations([1, 2, 5], 5), coin_permutations([1, 2, 5], 5)
(4, 9)

The four combinations are 5, 2+2+1, 2+1+1+1 and 1+1+1+1+1. The nine permutations count 2+2+1, 2+1+2 and 1+2+2 separately, and so on.

The reason: with coins in the outer loop, after processing coin k the array counts ways that use only the first k denominations. Each combination is built by deciding how many 1s, then how many 2s, then how many 5s, in that fixed order, so each multiset is generated exactly once. The intermediate arrays show it:

after coin 1   [1, 1, 1, 1, 1, 1]
after coin 2   [1, 1, 2, 2, 3, 3]
after coin 5   [1, 1, 2, 2, 3, 4]

With amounts in the outer loop, ways[a] sums over “which coin was last” with every denomination available at every step, so 1+2 and 2+1 arrive at 3 by different last coins and are counted separately. The array becomes [1, 1, 2, 3, 5, 9]. Climbing Stairs (1 or 2 steps at a time) is this permutation count with coins [1, 2], which is why it produces Fibonacci numbers. Combination Sum IV is the same count with arbitrary step sizes.

The forward inner loop in coin_combinations is correct here because coins are unlimited, the same reason it gives unbounded knapsack.

Subsequences and substrings

A subsequence keeps characters in order but may skip any of them (“ACE” is a subsequence of “ABCDE”). A substring is contiguous (“BCD”). The difference shows up directly in the recurrences: a subsequence DP can carry the best value forward past a mismatch, while a substring DP has to reset to zero.

Longest common subsequence

The longest common subsequence (LCS) of two strings is the longest sequence that is a subsequence of both. It’s the basis of diff, of version-control merge tools, and of sequence alignment.

• State: dp[i][j] is the LCS length of a[:i] and b[:j].
• Transition: if a[i-1] == b[j-1], that character can end the LCS, so dp[i][j] = dp[i-1][j-1] + 1. Otherwise at least one of the two last characters isn’t in the LCS, so dp[i][j] = max(dp[i-1][j], dp[i][j-1]).
• Base cases: dp[0][j] = dp[i][0] = 0.
• Answer: dp[m][n].

Taking the match when the last characters are equal is always safe. If an optimal LCS didn’t use that pair, it could be modified to end with it without getting shorter, which is the exchange argument that justifies not also checking the two max options in that case.

def lcs_table(a, b):
    m, n = len(a), len(b)
    dp = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if a[i - 1] == b[j - 1]:
                dp[i][j] = dp[i - 1][j - 1] + 1
            else:
                dp[i][j] = max(dp[i - 1][j], dp[i][j - 1])
    return dp


def lcs(a, b):
    dp = lcs_table(a, b)
    out, i, j = [], len(a), len(b)
    while i > 0 and j > 0:
        if a[i - 1] == b[j - 1]:
            out.append(a[i - 1])
            i, j = i - 1, j - 1
        elif dp[i - 1][j] >= dp[i][j - 1]:
            i -= 1
        else:
            j -= 1
    return "".join(reversed(out))

The table for a = "ABCBDAB" and b = "BDCABA":

         ""  B   D   C   A   B   A
    ""    0   0   0   0   0   0   0
    A     0   0   0   0   1   1   1
    B     0   1   1   1   1   2   2
    C     0   1   1   2   2   2   2
    B     0   1   1   2   2   3   3
    D     0   1   2   2   2   3   3
    A     0   1   2   2   3   3   4
    B     0   1   2   2   3   4   4

Reconstruction starts at the bottom-right 4 (row B, column A). The characters differ, and the cell above (row A, column A) also holds 4, so the walk moves up. There, A equals A: record “A” and move diagonally to (row D, column B), which holds 3. D and B differ; the cell above holds 3, move up to (row B, column B): a match, record “B”. Continuing the same way records “C” and then “B”. Reversed, the LCS is “BCBA”:

>>> lcs("ABCBDAB", "BDCABA")
'BCBA'

“BDAB” and “BCAB” are also common subsequences of length 4. Which one comes out depends on whether ties prefer moving up or left. Time and space are $O(mn)$. The length alone needs $O(\min(m, n))$ space with the two-row version shown earlier.

Longest common substring

For a contiguous match, the state has to say where the match ends. dp[i][j] is the length of the longest common suffix of a[:i] and b[:j], that is, the run of matching characters ending exactly at a[i-1] and b[j-1]. On a mismatch the run is broken, so the value is 0, not the max of the neighbors. The answer is the largest value anywhere in the table, not the bottom-right corner.

def longest_common_substring(a, b):
    best, end = 0, 0
    prev = [0] * (len(b) + 1)
    for i in range(1, len(a) + 1):
        cur = [0] * (len(b) + 1)
        for j in range(1, len(b) + 1):
            if a[i - 1] == b[j - 1]:
                cur[j] = prev[j - 1] + 1          # extend the run on the diagonal
                if cur[j] > best:
                    best, end = cur[j], i
        prev = cur
    return a[end - best:end]

>>> longest_common_substring("ABABC", "BABCA")
'BABC'

This is $O(mn)$ time and $O(n)$ space, and the substring itself is recovered from its end position and length, so the full table isn’t needed. For long strings, a suffix automaton or generalized suffix array solves the same problem in linear or near-linear time (chapter 5).

Longest increasing subsequence in O(n²)

The longest strictly increasing subsequence (LIS) of a sequence of numbers.

• State: length[i] is the length of the longest increasing subsequence that ends exactly at index i. The “ends exactly at” part is essential. “The LIS of the first i elements” can’t be extended, because extending requires knowing the last value.
• Transition: length[i] = 1 + max(length[j]) over j < i with nums[j] < nums[i], or 1 if there’s no such j.
• Answer: the maximum of length[i] over all i.

def lis_quadratic(nums):
    n = len(nums)
    length = [1] * n            # length[i] = LIS ending exactly at i
    parent = [-1] * n
    for i in range(n):
        for j in range(i):
            if nums[j] < nums[i] and length[j] + 1 > length[i]:
                length[i], parent[i] = length[j] + 1, j
    if not nums:
        return []
    i = max(range(n), key=length.__getitem__)
    out = []
    while i != -1:
        out.append(nums[i])
        i = parent[i]
    return out[::-1]

The parent array records which predecessor gave each length[i], so the subsequence is read back by following parents from the best ending position.

Longest increasing subsequence in O(n log n)

The faster method keeps a different invariant. tails[k] holds the smallest possible last value of any increasing subsequence of length k + 1 seen so far. A smaller tail is never worse than a larger one, because anything that extends the larger tail also extends the smaller.

tails is always strictly increasing, so it can be binary searched. For each new value x:

• If x is larger than every tail, it extends the longest subsequence: append it.
• Otherwise, find the first tail >= x and replace it with x. A subsequence of that length can now end on a smaller value.

The length of tails at the end is the LIS length. This is the card game patience sorting in disguise: each tail is the top card of a pile, a card goes on the leftmost pile whose top is not smaller, and the number of piles is the LIS length.

nums: 3 1 4 1 5 9 2 6 5 3 5 8 9 7 9

x=3  [3]
x=1  [1]                 1 replaces 3: a length-1 run can end lower
x=4  [1, 4]
x=1  [1, 4]              equal to tails[0]; strictly increasing, so replace with itself
x=5  [1, 4, 5]
x=9  [1, 4, 5, 9]
x=2  [1, 2, 5, 9]
x=6  [1, 2, 5, 6]
x=5  [1, 2, 5, 6]
x=3  [1, 2, 3, 6]
x=5  [1, 2, 3, 5]
x=8  [1, 2, 3, 5, 8]
x=9  [1, 2, 3, 5, 8, 9]
x=7  [1, 2, 3, 5, 7, 9]
x=9  [1, 2, 3, 5, 7, 9]

The final tails happens to be a valid subsequence here, but in general it is not: tails from different moments are mixed, and tails can contain a later small value in front of an earlier large one. To reconstruct a real subsequence, record for every element the index of the tail it was placed after at the moment it was placed:

from bisect import bisect_left

def lis(nums):
    tails = []        # tails[k] = smallest tail value of an increasing run of length k+1
    tail_idx = []     # index in nums of that tail
    parent = [-1] * len(nums)
    for i, x in enumerate(nums):
        k = bisect_left(tails, x)           # bisect_right would allow equal values
        if k == len(tails):
            tails.append(x)
            tail_idx.append(i)
        else:
            tails[k] = x
            tail_idx[k] = i
        parent[i] = tail_idx[k - 1] if k > 0 else -1
    out, i = [], tail_idx[-1] if tail_idx else -1
    while i != -1:
        out.append(nums[i])
        i = parent[i]
    return out[::-1]

>>> lis([3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9])
[1, 2, 3, 5, 7, 9]
>>> lis_quadratic([3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9])
[3, 4, 5, 6, 8, 9]

Both have length 6; they’re different optimal answers. bisect_left gives a strictly increasing subsequence. For non-decreasing (equal values allowed), use bisect_right.

Measured on random lists:

Four times the input took 16 times as long for the quadratic version and under five times as long for the patience method. A million elements in under a fifth of a second is out of reach for the quadratic version, which would need hours by extrapolation.

Classic problems built on LIS: Russian Doll Envelopes (sort by width ascending and, for equal widths, height descending, then take the LIS of heights; the descending tie-break stops two envelopes of the same width from nesting), building bridges between two banks, and the minimum number of deletions to make a sequence sorted (n - LIS).

Longest palindromic subsequence

The longest subsequence of a string that reads the same in both directions. This is the first interval DP in the chapter: the state is a substring s[i..j] identified by its two ends.

• State: dp[i][j] is the answer for s[i..j] (inclusive).
• Transition: if s[i] == s[j], both ends can wrap the best inner answer: dp[i+1][j-1] + 2. Otherwise one end is dropped: max(dp[i+1][j], dp[i][j-1]).
• Base case: dp[i][i] = 1.
• Order: dp[i][j] needs i + 1 and j - 1, so i runs downwards and j upwards (or, equivalently, by increasing interval length).

def longest_palindromic_subsequence(s):
    n = len(s)
    if n == 0:
        return 0
    dp = [[0] * n for _ in range(n)]
    for i in range(n - 1, -1, -1):
        dp[i][i] = 1
        for j in range(i + 1, n):
            if s[i] == s[j]:
                dp[i][j] = dp[i + 1][j - 1] + 2
            else:
                dp[i][j] = max(dp[i + 1][j], dp[i][j - 1])
    return dp[0][n - 1]

>>> longest_palindromic_subsequence("bbbab"), longest_palindromic_subsequence("character")
(4, 5)

“bbbb” and “carac” are the answers. When j = i + 1 and the characters match, dp[i+1][j-1] is dp[i+1][i], an empty interval whose value is the initial 0, so the result is 2 as it should be. Time and space are $O(n^2)$. The longest palindromic subsequence of s also equals the LCS of s and its reverse, which is a useful cross-check. The minimum number of insertions to make a string a palindrome is n minus this value.

String DP

Edit distance

The edit distance (Levenshtein distance, after Vladimir Levenshtein’s 1965 paper) between two strings is the minimum number of single-character insertions, deletions and substitutions that turn one into the other. Spell checkers rank suggestions by it, and fuzzy search in PostgreSQL’s fuzzystrmatch extension exposes it as levenshtein().

• State: dp[i][j] is the edit distance between a[:i] and b[:j].
• Transition: if a[i-1] == b[j-1], the last characters need no edit: dp[i-1][j-1]. Otherwise the last operation was one of three, each costing 1:
  • delete a[i-1]: dp[i-1][j] + 1
  • insert b[j-1]: dp[i][j-1] + 1
  • substitute a[i-1] with b[j-1]: dp[i-1][j-1] + 1
• Base cases: dp[i][0] = i (delete everything), dp[0][j] = j (insert everything).
• Answer: dp[m][n].

def edit_distance_table(a, b):
    m, n = len(a), len(b)
    dp = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(m + 1):
        dp[i][0] = i
    for j in range(n + 1):
        dp[0][j] = j
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if a[i - 1] == b[j - 1]:
                dp[i][j] = dp[i - 1][j - 1]
            else:
                dp[i][j] = 1 + min(dp[i - 1][j],       # delete a[i-1]
                                   dp[i][j - 1],       # insert b[j-1]
                                   dp[i - 1][j - 1])   # substitute
    return dp

The table for “kitten” to “sitting”:

         ""  s   i   t   t   i   n   g
    ""    0   1   2   3   4   5   6   7
    k     1   1   2   3   4   5   6   7
    i     2   2   1   2   3   4   5   6
    t     3   3   2   1   2   3   4   5
    t     4   4   3   2   1   2   3   4
    e     5   5   4   3   2   2   3   4
    n     6   6   5   4   3   3   2   3

The diagonal from (k, s) carries the 1 from the k-to-s substitution down through the matching “itt”. At (e, i) the value becomes 2 by substituting e with i. “n” matches “n” to give 2 at (n, n), and the final “g” has to be inserted, giving 3. Walking back from the corner and recording which neighbor produced each value gives the edit script:

def edit_script(a, b):
    dp = edit_distance_table(a, b)
    ops, i, j = [], len(a), len(b)
    while i > 0 or j > 0:
        if i > 0 and j > 0 and a[i - 1] == b[j - 1] and dp[i][j] == dp[i - 1][j - 1]:
            i, j = i - 1, j - 1
        elif i > 0 and j > 0 and dp[i][j] == dp[i - 1][j - 1] + 1:
            ops.append(f"substitute {a[i - 1]}->{b[j - 1]} at {i - 1}")
            i, j = i - 1, j - 1
        elif i > 0 and dp[i][j] == dp[i - 1][j] + 1:
            ops.append(f"delete {a[i - 1]} at {i - 1}")
            i -= 1
        else:
            ops.append(f"insert {b[j - 1]} at {i}")
            j -= 1
    return dp[-1][-1], ops[::-1]

>>> edit_script("kitten", "sitting")
(3, ['substitute k->s at 0', 'substitute e->i at 4', 'insert g at 6'])

The test harness applied each returned script to the source string and checked that it produced the target, in addition to checking the distance against a brute-force recursion.

Time $O(mn)$, space $O(mn)$, or $O(\min(m, n))$ for the distance alone. Variants change the costs or add operations. Allowing a swap of two adjacent characters gives the optimal string alignment distance, a restricted form of Damerau-Levenshtein. Weighted costs per character pair give the alignment scores used in bioinformatics (Needleman-Wunsch for global alignment, Smith-Waterman for local alignment, where cells are floored at zero in the same way that longest common substring resets to zero). LCS is the special case with insert and delete only: that distance is $m + n - 2 \cdot \text{LCS}$.

Word Break

Given a string and a dictionary, decide whether the string can be split into a sequence of dictionary words.

• State: ok[i] is True if the prefix s[:i] can be segmented.
• Transition: ok[i] is True if some j < i has ok[j] True and s[j:i] in the dictionary.
• Base case: ok[0] = True, the empty prefix.

def word_break(s, words):
    words = set(words)
    max_len = max(map(len, words), default=0)
    ok = [False] * (len(s) + 1)
    ok[0] = True                      # the empty prefix is always segmentable
    for i in range(1, len(s) + 1):
        for j in range(max(0, i - max_len), i):
            if ok[j] and s[j:i] in words:
                ok[i] = True
                break
    return ok[len(s)]

>>> word_break("applepenapple", ["apple", "pen"])
True
>>> word_break("catsandog", ["cats", "dog", "sand", "and", "cat"])
False

Limiting j to the last max_len positions bounds the work at $O(n \cdot L)$ slices of length at most $L$, so $O(n L^2)$ character work, instead of $O(n^3)$. A trie over the dictionary (chapter 11) avoids building slices: walk forward from each reachable j and mark every i where a word ends. The plain recursive solution without memoization is exponential on inputs like "aaaa...ab" with dictionary ["a", "aa", "aaa"], which is the standard demonstration of why the cache matters. Word Break II (list every segmentation) is a backtracking problem (chapter 15) whose output can be exponential in size, though the ok array is still useful for pruning dead branches.

Palindrome partitioning: minimum cuts

Cut a string into pieces that are all palindromes, using as few cuts as possible. Two tables work together:

• is_pal[i][j] says whether s[i..j] is a palindrome. It’s an interval DP: the ends match and the inside is a palindrome.
• cuts[j] is the fewest cuts for s[:j+1]. If the whole prefix is a palindrome it’s 0. Otherwise the last piece is some palindrome s[i..j], and the cost is cuts[i-1] + 1.

def min_palindrome_cuts(s):
    n = len(s)
    is_pal = [[False] * n for _ in range(n)]
    for i in range(n - 1, -1, -1):
        for j in range(i, n):
            is_pal[i][j] = s[i] == s[j] and (j - i < 2 or is_pal[i + 1][j - 1])
    cuts = [0] * n                    # cuts[j] = min cuts for s[:j+1]
    for j in range(n):
        if is_pal[0][j]:
            continue
        cuts[j] = min(cuts[i - 1] + 1 for i in range(1, j + 1) if is_pal[i][j])
    return cuts[-1] if n else 0

>>> min_palindrome_cuts("aab"), min_palindrome_cuts("banana"), min_palindrome_cuts("abcde")
(1, 1, 4)

“aa|b”, “b|anana”, and four cuts between five distinct letters. Both tables are $O(n^2)$ in time and the palindrome table is $O(n^2)$ in space. A version that expands around each center instead of storing is_pal brings space down to $O(n)$. The minimum-cut version is DP; listing every palindrome partition is backtracking.

Wildcard matching

The pattern has ?, which matches any single character, and *, which matches any sequence including the empty one, as in shell globs.

• State: dp[i][j] is True if s[:i] matches p[:j].
• Transition: for a literal or ?, the last characters must match and dp[i-1][j-1] must hold. For *, either it matches nothing (dp[i][j-1]) or it absorbs s[i-1] and stays available for more (dp[i-1][j]).
• Base cases: dp[0][0] = True. An empty string matches a pattern prefix only if that prefix is all stars.

def wildcard_match(s, p):
    m, n = len(s), len(p)
    dp = [[False] * (n + 1) for _ in range(m + 1)]
    dp[0][0] = True
    for j in range(1, n + 1):
        dp[0][j] = dp[0][j - 1] and p[j - 1] == "*"
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if p[j - 1] == "*":
                dp[i][j] = dp[i][j - 1] or dp[i - 1][j]   # * matches empty, or one more char
            elif p[j - 1] == "?" or p[j - 1] == s[i - 1]:
                dp[i][j] = dp[i - 1][j - 1]
    return dp[m][n]

>>> wildcard_match("adceb", "*a*b"), wildcard_match("acdcb", "a*c?b")
(True, False)

The dp[i-1][j] term for * is the key step. It lets one star consume any number of characters in $O(1)$ per cell instead of looping over every possible length, keeping the total at $O(mn)$. Results were checked against fnmatch.fnmatchcase on random strings and patterns.

Regular expression matching

Here . matches any one character and x* means zero or more copies of the preceding element x. The star now binds to the character before it, which changes the transitions:

• p[j-1] is *: either use zero copies of p[j-2], dropping both pattern characters (dp[i][j-2]), or, if p[j-2] matches s[i-1], let x* absorb that character and stay in place (dp[i-1][j]).
• Otherwise: p[j-1] must be . or equal s[i-1], and dp[i-1][j-1] must hold.

def regex_match(s, p):
    m, n = len(s), len(p)
    dp = [[False] * (n + 1) for _ in range(m + 1)]
    dp[0][0] = True
    for j in range(2, n + 1):
        dp[0][j] = p[j - 1] == "*" and dp[0][j - 2]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if p[j - 1] == "*":
                zero = dp[i][j - 2]                                   # drop "x*"
                more = p[j - 2] in (s[i - 1], ".") and dp[i - 1][j]  # x* eats s[i-1]
                dp[i][j] = zero or more
            elif p[j - 1] in (s[i - 1], "."):
                dp[i][j] = dp[i - 1][j - 1]
    return dp[m][n]

>>> regex_match("aab", "c*a*b"), regex_match("mississippi", "mis*is*p*."), regex_match("ab", ".*")
(True, False, True)

The code assumes a well-formed pattern in which * never appears first. It was checked against re.fullmatch on random strings and patterns built from a, b, . and starred elements. A backtracking matcher for the same syntax can take exponential time on patterns like a*a*a*a*c against a long run of as. The DP is $O(mn)$ regardless, which is the same guarantee that automaton-based regex engines provide for full regular expressions.

Grid DP

In grid problems the state is a cell, and movement restrictions (usually only right and down) make the dependency graph acyclic, so cells can be filled row by row.

Unique paths

Count the paths from the top-left to the bottom-right of an r by c grid moving only right or down. Every cell is reached from the cell above or the cell to the left, so paths[r][c] = paths[r-1][c] + paths[r][c-1], with 1 along the top row and left column. For a 3 by 7 grid:

1   1   1   1   1   1   1
1   2   3   4   5   6   7
1   3   6  10  15  21  28

The table is Pascal’s triangle turned on its side. Without obstacles there’s a closed form: a path is a sequence of $r - 1$ downs and $c - 1$ rights, so the count is $\binom{r + c - 2}{r - 1}$, and $\binom{8}{2} = 28$. The DP becomes necessary once obstacles appear:

def unique_paths_obstacles(grid):
    """grid[r][c] == 1 means blocked. Moves are right and down only."""
    cols = len(grid[0])
    row = [0] * cols
    row[0] = 1 if grid[0][0] == 0 else 0
    for r in range(len(grid)):
        for c in range(cols):
            if grid[r][c] == 1:
                row[c] = 0
            elif c > 0:
                row[c] += row[c - 1]      # row[c] still holds the value from above
    return row[-1]

>>> unique_paths_obstacles([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
2

This is the one-row rolling array in its most compact form. Before row[c] is updated it holds the count from the cell above, and row[c-1] has already been updated to the current row’s value from the left, so a left-to-right sweep gives exactly the needed sum. Time $O(rc)$, space $O(c)$. A blocked start or end correctly yields 0. The tests compared it with explicit enumeration of every path on random grids, and with the binomial formula on grids without obstacles.

Minimum path sum

Each cell has a non-negative cost. Find the cheapest path from top-left to bottom-right moving right or down. Replace the sum with a minimum and add the cell’s own cost:

def min_path_sum(grid):
    rows, cols = len(grid), len(grid[0])
    dp = [[0] * cols for _ in range(rows)]
    for r in range(rows):
        for c in range(cols):
            if r == 0 and c == 0:
                dp[r][c] = grid[r][c]
            elif r == 0:
                dp[r][c] = dp[r][c - 1] + grid[r][c]
            elif c == 0:
                dp[r][c] = dp[r - 1][c] + grid[r][c]
            else:
                dp[r][c] = min(dp[r - 1][c], dp[r][c - 1]) + grid[r][c]
    path, r, c = [(rows - 1, cols - 1)], rows - 1, cols - 1
    while (r, c) != (0, 0):
        if c == 0 or (r > 0 and dp[r - 1][c] <= dp[r][c - 1]):
            r -= 1
        else:
            c -= 1
        path.append((r, c))
    return dp[-1][-1], path[::-1]

>>> min_path_sum([[1, 3, 1], [1, 5, 1], [4, 2, 1]])
(7, [(0, 0), (0, 1), (0, 2), (1, 2), (2, 2)])

The path follows 1, 3, 1, 1, 1. If moves in all four directions are allowed, the dependency graph has cycles and the problem becomes a shortest-path problem for Dijkstra’s algorithm (chapter 13). DP requires an acyclic order of states.

Dungeon Game, and why it runs backwards

A knight starts at the top-left of a grid and must reach a princess at the bottom-right, moving right or down. Each cell adds to or subtracts from his health, and if health ever drops to 0 he dies. Find the minimum starting health.

The natural forward DP fails. A forward state “best health on arrival at this cell” has to trade off two quantities that don’t combine: the current health and the lowest health along the way so far. A path that has kept health high and a path that has collected a lot of health but dipped dangerously early are incomparable, so keeping only one number per cell loses information and gives wrong answers.

Running backwards removes the problem, because the future fully determines what’s needed. Define need[r][c] as the minimum health required on entering cell (r, c) to survive from there to the end. Leaving the cell, the knight needs after = min(need[below], need[right]). Entering with health h, he has h + grid[r][c] after the cell, which must be at least after, and h itself must be at least 1:

$\text{need}[r][c] = \max\big(1,\ \min(\text{need}[r+1][c],\ \text{need}[r][c+1]) - \text{grid}[r][c]\big)$

def dungeon_min_health(grid):
    rows, cols = len(grid), len(grid[0])
    need = [[inf] * (cols + 1) for _ in range(rows + 1)]
    need[rows][cols - 1] = need[rows - 1][cols] = 1      # health 1 after the last room
    for r in range(rows - 1, -1, -1):
        for c in range(cols - 1, -1, -1):
            after = min(need[r + 1][c], need[r][c + 1])
            need[r][c] = max(1, after - grid[r][c])
    return need[0][0]

For the standard example:

grid                    need (computed from the bottom-right)
 -2   -3    3            7    5    2
 -5  -10    1            6   11    5
 10   30   -5            1    1    6

>>> dungeon_min_health([[-2, -3, 3], [-5, -10, 1], [10, 30, -5]])
7

Starting with 7, the path right, right, down, down goes 7 → 5 → 2 → 5 → 6 → 1. The max(1, ...) clamp is what makes backwards work: a big health potion later can’t retroactively save the knight from dying earlier, and the clamp encodes that a cell can never require less than 1 on entry, no matter how much it gives. The test compared the result with the true answer computed for every path as one minus the lowest prefix sum.

The general lesson: if the constraint is about what must hold along the way and the objective depends on the future, try defining the state as “what is needed from here to the end”.

Maximum sum rectangle

Find the sub-rectangle of a matrix with the largest sum. The one-dimensional version, the maximum sum subarray, is solved by Kadane’s algorithm, itself a tiny DP: cur is the best sum of a subarray ending at the current position, and it either extends the previous run or restarts at the current element.

def kadane(arr):
    best = cur = arr[0]
    for x in arr[1:]:
        cur = max(x, cur + x)
        best = max(best, cur)
    return best

For 2D, fix a top row and a bottom row. Summing each column between them collapses the rectangle to a 1D array, whose best subarray is the best rectangle with exactly those top and bottom rows. Extending bottom by one row only adds that row to the column sums, so each pair of rows costs $O(\text{cols})$:

def max_sum_rectangle(matrix):
    rows, cols = len(matrix), len(matrix[0])
    best = -inf
    for top in range(rows):
        col_sums = [0] * cols
        for bottom in range(top, rows):
            for c in range(cols):
                col_sums[c] += matrix[bottom][c]   # sums of rows top..bottom per column
            best = max(best, kadane(col_sums))
    return best

>>> max_sum_rectangle([[1, 2, -1, -4, -20], [-8, -3, 4, 2, 1], [3, 8, 10, 1, 3], [-4, -1, 1, 7, -6]])
29

The answer is the 3-row, 3-column block starting at row 1, column 1 (-3 4 2 / 8 10 1 / -1 1 7). Time is $O(R^2 C)$, so put the smaller dimension in the outer pair of loops. Brute force over all rectangles with direct summation is $O(R^3 C^3)$, and with 2D prefix sums $O(R^2 C^2)$. Kadane’s version, initialized with the first element rather than 0, returns the largest single element when every number is negative, which is correct when the rectangle must be non-empty.

Other patterns

Rod cutting

A rod of length n can be cut into integer-length pieces, and a piece of length k sells for prices[k]. Find the most profitable way to cut it. The last piece (or equivalently the first) has some length k, and the rest is a smaller rod:

$\text{best}[n] = \max_{1 \le k \le n} \big(\text{prices}[k] + \text{best}[n - k]\big)$

It’s an unbounded knapsack in which the capacity is the rod length and every piece length is an item. Storing the k that achieved each maximum gives the cuts:

def rod_cutting(prices, n):
    """prices[k] is the price of a piece of length k (prices[0] unused)."""
    best = [0] + [-inf] * n             # -inf so the first candidate always wins
    first_cut = [0] * (n + 1)
    for length in range(1, n + 1):
        for k in range(1, min(length, len(prices) - 1) + 1):
            if prices[k] + best[length - k] > best[length]:
                best[length] = prices[k] + best[length - k]
                first_cut[length] = k
    pieces, rest = [], n
    while rest > 0:
        pieces.append(first_cut[rest])
        rest -= first_cut[rest]
    return best[n], pieces

With the prices from Introduction to Algorithms (Cormen, Leiserson, Rivest and Stein), [0, 1, 5, 8, 9, 10, 17, 17, 20, 24, 30]:

>>> prices = [0, 1, 5, 8, 9, 10, 17, 17, 20, 24, 30]
>>> for n in (4, 7, 8, 10):
...     print(n, rod_cutting(prices, n))
...
4 (10, [2, 2])
7 (18, [1, 6])
8 (22, [2, 6])
10 (30, [10])

Initializing best to minus infinity rather than 0 matters for reconstruction. With 0 and a price list of all zeros, no candidate is ever strictly better, first_cut stays 0, and the reconstruction loop never terminates. The first version of this code had exactly that bug, and the random tests found it.

Time $O(n^2)$, space $O(n)$. The naive recursion that tries every first cut is $O(2^{n-1})$, one call per composition of n. The code assumes every piece length up to n has a price; the min(...) bound allows a shorter price list, and the test cases included rods longer than the price list.

Matrix chain multiplication (interval DP)

Multiplying a $p \times q$ matrix by a $q \times r$ matrix takes $pqr$ scalar multiplications. Matrix multiplication is associative, so a chain $A_1 A_2 \cdots A_n$ can be parenthesized in any way with the same result, but the costs differ enormously. With shapes $10 \times 30$, $30 \times 5$ and $5 \times 60$:

• $(A_1 A_2) A_3$ costs $10 \cdot 30 \cdot 5 + 10 \cdot 5 \cdot 60 = 1{,}500 + 3{,}000 = 4{,}500$.
• $A_1 (A_2 A_3)$ costs $30 \cdot 5 \cdot 60 + 10 \cdot 30 \cdot 60 = 9{,}000 + 18{,}000 = 27{,}000$.

The number of parenthesizations is a Catalan number, which grows like $4^n / n^{3/2}$, so trying them all is hopeless beyond a dozen matrices. This is the canonical interval DP: the state is a contiguous range, and the transition tries every split point inside it.

• State: cost[i][j] is the cheapest way to multiply matrices i through j, where matrix i has shape dims[i] × dims[i+1].
• Transition: the last multiplication joins the product of i..k with the product of k+1..j. Their shapes are dims[i] × dims[k+1] and dims[k+1] × dims[j+1], so $\text{cost}[i][j] = \min_{i \le k < j} \big(\text{cost}[i][k] + \text{cost}[k+1][j] + d_i\, d_{k+1}\, d_{j+1}\big)$
• Base case: cost[i][i] = 0.
• Order: by increasing interval length, so both halves of every split are ready.

def matrix_chain(dims):
    """Matrix i has shape dims[i] x dims[i+1]. Returns (cost, parenthesization)."""
    n = len(dims) - 1
    cost = [[0] * n for _ in range(n)]
    split = [[0] * n for _ in range(n)]
    for length in range(2, n + 1):              # interval length, short to long
        for i in range(n - length + 1):
            j = i + length - 1
            cost[i][j] = inf
            for k in range(i, j):
                c = cost[i][k] + cost[k + 1][j] + dims[i] * dims[k + 1] * dims[j + 1]
                if c < cost[i][j]:
                    cost[i][j], split[i][j] = c, k

    def show(i, j):
        if i == j:
            return f"A{i + 1}"
        k = split[i][j]
        return f"({show(i, k)} {show(k + 1, j)})"

    return cost[0][n - 1], show(0, n - 1)

>>> matrix_chain([10, 30, 5, 60])
(4500, '((A1 A2) A3)')
>>> matrix_chain([30, 35, 15, 5, 10, 20, 25])
(15125, '((A1 (A2 A3)) ((A4 A5) A6))')

The second is the six-matrix example from Cormen et al. Time $O(n^3)$ ($O(n^2)$ intervals, $O(n)$ splits each), space $O(n^2)$. The tests compared the cost against explicit enumeration of every bracketing, and also re-evaluated the returned parenthesization string to confirm it achieves the claimed cost.

Iterating by interval length is the pattern to remember. Iterating i and j in plain row order would read cost[k+1][j] before it’s computed. The same skeleton solves Burst Balloons, optimal binary search tree construction, polygon triangulation, and the minimum cost to merge stones. Some interval DPs with a monotone optimal split point can be sped up to $O(n^2)$ with Knuth’s optimization (it applies to optimal binary search trees), but matrix chain doesn’t satisfy its conditions in general. Hu and Shing published an $O(n \log n)$ algorithm for matrix chain specifically, which is rarely implemented.

Egg drop

With k identical eggs and a building of n floors, there is an unknown critical floor f between 0 and n: an egg dropped from above f breaks, at or below f it doesn’t. A broken egg can’t be reused. Find the minimum number of drops that determines f in the worst case.

The classic DP. Let dp[e][f] be the worst-case drops for e eggs and f candidate floors. Dropping from floor x either breaks the egg (search the x - 1 floors below with e - 1 eggs) or doesn’t (search the f - x floors above with e eggs). The adversary picks the worse outcome and the solver picks the best x:

$\text{dp}[e][f] = 1 + \min_{1 \le x \le f} \max\big(\text{dp}[e-1][x-1],\ \text{dp}[e][f-x]\big)$

with dp[1][f] = f (one egg forces a linear scan from the bottom) and dp[e][0] = 0.

def egg_drop_classic(eggs, floors):
    # dp[e][f] = worst-case drops to find the critical floor among f floors with e eggs
    dp = [[0] * (floors + 1) for _ in range(eggs + 1)]
    for f in range(1, floors + 1):
        dp[1][f] = f
    for e in range(2, eggs + 1):
        for f in range(1, floors + 1):
            dp[e][f] = min(1 + max(dp[e - 1][x - 1], dp[e][f - x])
                           for x in range(1, f + 1))
    return dp[eggs][floors]

This is $O(k n^2)$ time. Because the first term of the max increases with x and the second decreases, the best x can be binary searched, giving $O(k n \log n)$.

The moves-based formulation. Flip the question: with m drops and e eggs, how many floors can be fully resolved? Call it covered[m][e]. The first drop is made from some floor. If the egg breaks, the m - 1 remaining drops and e - 1 eggs can resolve covered[m-1][e-1] floors below it. If it survives, the same drops with e eggs resolve covered[m-1][e] floors above. Counting the floor dropped from:

$\text{covered}[m][e] = \text{covered}[m-1][e-1] + \text{covered}[m-1][e] + 1$

The answer is the smallest m with covered[m][k] >= n. A single array updated from high e to low (the same backwards trick as 0/1 knapsack, so that covered[e-1] still holds the previous move’s value) is enough:

def egg_drop_moves(eggs, floors):
    # covered[e] = most floors that can be resolved with e eggs in the current number of moves
    covered = [0] * (eggs + 1)
    moves = 0
    while covered[eggs] < floors:
        moves += 1
        for e in range(eggs, 0, -1):
            covered[e] = covered[e - 1] + covered[e] + 1
    return moves

The first few moves for 1, 2 and 3 eggs:

moves   1 egg   2 eggs   3 eggs
  1       1       1        1
  2       2       3        3
  3       3       6        7
  4       4      10       14
  5       5      15       25
  6       6      21       41
  7       7      28       63

With one egg the coverage grows linearly, with two eggs it’s the triangular numbers $m(m+1)/2$, and with unlimited eggs it reaches $2^m - 1$, which is binary search. In general $\text{covered}[m][e] = \sum_{i=1}^{e} \binom{m}{i}$.

>>> egg_drop_classic(2, 100), egg_drop_moves(2, 100), egg_drop_moves(3, 100)
(14, 14, 9)
>>> egg_drop_moves(2, 10**9), egg_drop_moves(20, 10**18)
(44721, 75)

The moves-based version runs in $O(k \cdot m)$ where $m$ is the answer, which is at most n and usually far smaller, and it handles $10^{18}$ floors instantly. Both formulations agreed with each other and with a direct memoized minimax on every combination of up to 5 eggs and 80 floors. Measured:

The lesson generalizes. When the answer is small but one of the inputs is huge, swap them: make the answer part of the state and compute the largest input that answer can handle.

DP on trees

On a rooted tree, a node’s subtree is a natural subproblem. A post-order traversal computes each node’s value from its children’s values, so the evaluation order comes for free. The common move is to return more than one value per node, one for each relevant situation of the node itself.

House Robber III. The houses form a binary tree, and robbing a node and its parent triggers the alarm. Each node returns a pair: the best total in its subtree if it is robbed, and if it isn’t.

class Node:
    def __init__(self, val, left=None, right=None):
        self.val, self.left, self.right = val, left, right


def rob_tree(root):
    def go(node):
        if node is None:
            return 0, 0                       # (best if robbed, best if not)
        l_rob, l_skip = go(node.left)
        r_rob, r_skip = go(node.right)
        robbed = node.val + l_skip + r_skip
        skipped = max(l_rob, l_skip) + max(r_rob, r_skip)
        return robbed, skipped
    return max(go(root))

>>> root = Node(3, Node(2, None, Node(3)), Node(3, None, Node(1)))
>>> rob_tree(root)
7

If the node is robbed, its children must not be. If it isn’t, each child independently takes whichever of its two values is larger. This is maximum-weight independent set, NP-hard on general graphs but $O(n)$ on trees because subtrees don’t interact. The recursion depth equals the tree height, so a degenerate (path-like) tree of more than about a thousand nodes needs either a higher recursion limit or an iterative post-order traversal.

Tree diameter. The longest path between any two nodes. Every path has a single highest node, where it bends, so for each node u take the two longest downward paths into different children and add them:

def tree_diameter(n, edges):
    """Longest path (in edges) in an unrooted tree with nodes 0..n-1."""
    adj = [[] for _ in range(n)]
    for u, v in edges:
        adj[u].append(v)
        adj[v].append(u)
    best = 0
    depth = [0] * n                     # longest downward path from each node
    order, parent = [0], [-1] * n
    for u in order:                     # BFS order; children come after parents
        for v in adj[u]:
            if v != parent[u]:
                parent[v] = u
                order.append(v)
    for u in reversed(order):           # children before parents
        top1 = top2 = 0
        for v in adj[u]:
            if v != parent[u]:
                d = depth[v] + 1
                if d > top1:
                    top1, top2 = d, top1
                elif d > top2:
                    top2 = d
        depth[u] = top1
        best = max(best, top1 + top2)   # path bending at u
    return best

Processing nodes in reverse BFS order visits every child before its parent, which gives a post-order without recursion. The two-BFS method (BFS from any node to the farthest node x, then BFS from x; the largest distance found is the diameter) is equally fast for unweighted trees, but the DP version extends directly to weighted edges, including negative ones, where the two-BFS argument fails.

Rerooting. Some tree questions ask for an answer for every node as root, for example the sum of distances from each node to all others (LeetCode’s Sum of Distances in Tree). Running a separate traversal from every node costs $O(n^2)$. Rerooting does it in two passes, $O(n)$ total:

1. A post-order pass from an arbitrary root computes, for every node, its subtree size and down[u], the sum of distances from u to the nodes in its own subtree.
2. A pre-order pass moves the root from a parent p to its child u. The size[u] nodes in u’s subtree each get one step closer, and the other n - size[u] nodes each get one step farther, so answer[u] = answer[p] - size[u] + (n - size[u]).

def sum_of_distances(n, edges):
    """For every node, the sum of distances to all other nodes. O(n) by rerooting."""
    adj = [[] for _ in range(n)]
    for u, v in edges:
        adj[u].append(v)
        adj[v].append(u)
    order, parent = [0], [-1] * n
    for u in order:
        for v in adj[u]:
            if v != parent[u]:
                parent[v] = u
                order.append(v)
    size = [1] * n
    down = [0] * n                      # sum of distances from u into its subtree
    for u in reversed(order):
        if parent[u] != -1:
            p = parent[u]
            size[p] += size[u]
            down[p] += down[u] + size[u]
    answer = [0] * n
    answer[0] = down[0]
    for u in order[1:]:                 # parents before children
        # moving the root from parent to u: size[u] nodes get 1 closer, the rest 1 farther
        answer[u] = answer[parent[u]] - size[u] + (n - size[u])
    return answer

>>> edges = [(0, 1), (0, 2), (2, 3), (2, 4), (2, 5)]
>>> tree_diameter(6, edges), sum_of_distances(6, edges)
(3, [8, 12, 6, 10, 10, 10])

Both were compared with a BFS from every node on random trees of up to 30 nodes. The rerooting pattern applies whenever the answer at a child can be derived from the answer at its parent by removing the child’s own contribution and adding the “rest of the tree” contribution. For operations without an inverse (like max), store prefix and suffix maxima over each node’s children so a child’s contribution can be excluded without subtraction.

DP with bitmasks

When the state must record which of a small set of items have been used, a bitmask encodes the subset as an integer: bit i set means item i is in it. With n items there are $2^n$ masks, which is practical up to about n = 20 in compiled languages, and somewhat less in Python.

Traveling salesman (Held-Karp). Find the shortest tour that visits every city exactly once and returns to the start. Trying all orders costs $O(n!)$. The Held-Karp algorithm (Bellman, and independently Held and Karp, both 1962) notices that the cost of finishing a tour depends only on the set of cities already visited and the current city, not on the order in which they were visited.

• State: dp[mask][last] is the shortest path that starts at city 0, visits exactly the cities in mask, and ends at last.
• Transition: extend to any unvisited nxt: dp[mask | 1 << nxt][nxt] = min(..., dp[mask][last] + dist[last][nxt]).
• Base case: dp[{0}][0] = 0.
• Order: increasing mask, since adding a bit always produces a larger integer.
• Answer: min(dp[all][j] + dist[j][0]).

def held_karp(dist):
    """Shortest tour visiting every city once, starting and ending at city 0."""
    n = len(dist)
    full = 1 << n
    dp = [[inf] * n for _ in range(full)]
    parent = [[-1] * n for _ in range(full)]
    dp[1][0] = 0                        # visited {0}, standing at 0
    for mask in range(1, full, 2):      # only masks that contain city 0
        for last in range(n):
            cost = dp[mask][last]
            if cost == inf:
                continue
            for nxt in range(n):
                if mask & (1 << nxt):
                    continue
                new_mask = mask | (1 << nxt)
                if cost + dist[last][nxt] < dp[new_mask][nxt]:
                    dp[new_mask][nxt] = cost + dist[last][nxt]
                    parent[new_mask][nxt] = last
    mask = full - 1
    best, last = min((dp[mask][j] + dist[j][0], j) for j in range(n))
    tour = []
    while last != -1:
        tour.append(last)
        mask, last = mask ^ (1 << last), parent[mask][last]
    return best, tour[::-1] + [0]

>>> d = [[0, 10, 15, 20], [10, 0, 35, 25], [15, 35, 0, 30], [20, 25, 30, 0]]
>>> held_karp(d)
(80, [0, 2, 3, 1, 0])

Time $O(2^n n^2)$, space $O(2^n n)$. It doesn’t assume symmetric distances or the triangle inequality. Compared with trying every permutation:

At 10 cities Held-Karp is already far ahead, and the gap grows with every city, since $n!$ outruns $2^n n^2$ quickly. At 14 cities there are $13! \approx 6.2 \times 10^9$ permutations; extrapolating from the 10-city run, about 0.5 µs each, that’s on the order of an hour in Python, against 51 ms for Held-Karp. Memory becomes the limit next: at n = 20 the table has about 21 million entries.

Assignment-style subset DP. Assign n workers to n jobs, one each, minimizing total cost. If workers are assigned in a fixed order, the number of bits set in the mask tells whose turn it is, so the state is just the mask:

def min_assignment(cost):
    """Assign n workers to n jobs, one each, minimizing total cost. O(2^n * n)."""
    n = len(cost)
    dp = [inf] * (1 << n)               # dp[mask] = best cost giving jobs in mask
    dp[0] = 0                           # to the first popcount(mask) workers
    for mask in range(1 << n):
        if dp[mask] == inf:
            continue
        worker = mask.bit_count()
        if worker == n:
            continue
        for job in range(n):
            if not mask & (1 << job):
                nxt = mask | (1 << job)
                dp[nxt] = min(dp[nxt], dp[mask] + cost[worker][job])
    return dp[-1]

>>> min_assignment([[9, 2, 7, 8], [6, 4, 3, 7], [5, 8, 1, 8], [7, 6, 9, 4]])
13

That’s worker 0 to job 1 (2), worker 1 to job 0 (6), worker 2 to job 2 (1) and worker 3 to job 3 (4). Time $O(2^n n)$ versus $O(n!)$ for permutations. For this particular problem the Hungarian algorithm is polynomial, $O(n^3)$, and is the right choice for large n. The bitmask DP earns its place when the cost isn’t a simple sum of independent pairs, for instance when it depends on which jobs were already taken, which the Hungarian algorithm can’t model. Related problems: Partition to K Equal Sum Subsets, minimum number of work sessions, Smallest Sufficient Team, and the “iterate over submasks” technique (sub = (sub - 1) & mask) whose total cost over all masks is $O(3^n)$.

Digit DP

Digit DP counts the integers in a range [lo, hi] with some property of their digits, when the range is far too large to scan. The standard reduction is count(hi) - count(lo - 1), where count(N) counts qualifying integers in [0, N]. Then count(N) builds numbers digit by digit from the most significant end, carrying:

• pos: the current digit position.
• tight: whether every digit so far equals N’s digit at that position. While tight, the next digit can go only up to N’s digit; once a smaller digit has been placed, any later digit up to 9 is allowed and the number is guaranteed to stay below N.
• The property’s own state, here the digit sum modulo k.

The example counts integers whose digit sum is divisible by k:

def count_digit_sum_divisible(limit, k):
    """How many integers x with 0 <= x <= limit have digit sum divisible by k."""
    digits = list(map(int, str(limit)))

    @cache
    def go(pos, rem, tight):
        if pos == len(digits):
            return 1 if rem == 0 else 0
        top = digits[pos] if tight else 9
        return sum(go(pos + 1, (rem + d) % k, tight and d == top)
                   for d in range(top + 1))

    return go(0, 0, True)


def count_in_range(lo, hi, k):
    below = count_digit_sum_divisible(lo - 1, k) if lo > 0 else 0
    return count_digit_sum_divisible(hi, k) - below

>>> count_in_range(1, 100, 7), count_in_range(1, 10**18, 7)
(13, 142857142856594041)

Shorter numbers are handled implicitly as numbers with leading zeros, which don’t change a digit sum. That’s convenient here but must be handled with an extra started flag for properties where leading zeros matter, such as “no two adjacent digits are equal” or “contains the digit 0”. The number of states is $(\text{digits}) \times k \times 2$ and each tries 10 digits, so the $10^{18}$ range takes about $19 \times 7 \times 2 \times 10$ steps. The function was checked against a direct scan of every integer for random ranges up to 6,000 and random k up to 12. Note that the @cache lives inside the function, so each call gets a fresh cache tied to its own digits. A module-level cache keyed only on (pos, rem, tight) would return answers for the wrong limit.

Catalan numbers

The Catalan numbers 1, 1, 2, 5, 14, 42, 132, … satisfy the recurrence

$C_0 = 1, \qquad C_{n} = \sum_{i=0}^{n-1} C_i\, C_{n-1-i}$

and the closed form $C_n = \frac{1}{n+1}\binom{2n}{n}$. The recurrence is itself a DP. For binary search trees on keys $1..n$: choose the root $i + 1$, which leaves $i$ keys for the left subtree and $n - 1 - i$ for the right, and the two sides are independent.

def catalan(n):
    c = [1] + [0] * n
    for m in range(1, n + 1):
        c[m] = sum(c[i] * c[m - 1 - i] for i in range(m))   # root splits the rest
    return c

>>> catalan(10)
[1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796]

The DP takes $O(n^2)$ and was checked against math.comb(2n, n) // (n + 1) and against a brute-force count of balanced parenthesis strings. $C_n$ counts, among many other things:

• structurally distinct binary search trees with $n$ keys (Unique Binary Search Trees),
• strings of $n$ pairs of balanced parentheses,
• full binary trees with $n + 1$ leaves, and therefore the ways to parenthesize a product of $n + 1$ factors (the matrix chain search space),
• triangulations of a convex polygon with $n + 2$ sides,
• monotone lattice paths across an $n \times n$ grid that never cross the diagonal (Dyck paths),
• sequences of $n$ pushes and $n$ pops that are valid for a stack.

When a counting problem splits into “a first piece, then an independent left part and right part whose sizes add up to $n - 1$”, Catalan numbers are usually the answer. For a single value, the closed form is faster. The DP pattern matters because it adapts to variants the closed form doesn’t cover, such as counting trees with a height limit.