Algorithm design: divide and conquer, greedy, dynamic programming, and flow

Why this subject matters now

The complaint that classical algorithms stopped mattering once machines got fast is exactly backwards, and the last few years have made the argument for the subject sharper than it has been in a while. Three things changed. First, the constants moved but the exponents did not. A modern server does on the order of \(10^{10}\) useful operations per second per core, which means an \(O(n^2)\) routine on a ten-million-element input is a machine-week and an \(O(n \log n)\) routine on the same input is a fifth of a second. No amount of hardware progress closes a gap that is a function of \(n\), and inputs have grown faster than clocks for two decades. Second, the interesting resource stopped being instructions and became memory traffic and parallelism, which changed which algorithms win without changing how they are analyzed. The reason a fused attention kernel beats a materialized one is an asymptotic argument about bytes moved, not FLOPs, and the reason scan-based and tree-based formulations dominate on accelerators is a depth-versus-work analysis of the same recurrences solved below. Third, and least appreciated, the practical meaning of NP-hardness inverted. In 1990 an NP-hard model was a dead end. Today the combination of conflict-driven clause learning in SAT solvers, branch-and-cut in mixed-integer programming, and specialized solvers such as OR-Tools' CP-SAT means an engineer who can express a scheduling or routing or assignment problem as an integer program or a SAT instance routinely gets exact answers to instances with hundreds of thousands of variables. Modeling ability, meaning the ability to see one problem as another, is now worth more than implementation ability, and modeling ability is precisely what proofs of reduction teach.

There is a fourth change that matters for anybody working near machine learning. Program synthesis by language models has made writing a correct implementation of a known algorithm nearly free, while leaving completely untouched the two skills that decide whether the resulting system works, choosing the right abstraction (is this a flow problem, a matroid, a DP over subsets, or genuinely hard) and proving the bound (does this terminate, is the greedy choice safe, what is the state space). The valuable half of the subject is the half that is hardest to automate, which is why this page spends its length on derivations and reductions and gives implementations only where an implementation settles something. Median-of-medians appears because the pivot argument is subtle, the FFT because the recursion is easy to get subtly wrong, Dinic's algorithm because it is the workhorse behind every flow reduction, and the dynamic programs because the recurrence is the artifact and the code is a transcription of it.

Core theory

Analysis foundations

The RAM model and what it assumes

Running time has to be counted against a machine model, and the standard one is the word RAM. It postulates a processor with a constant number of registers and an unbounded array of memory cells, each holding a \(w\)-bit word. It fixes an instruction set of arithmetic (add, subtract, multiply, divide, compare), bitwise operations, and indexed load and store, each costing one unit of time, and it adds the crucial assumption \( w = \Theta(\log n) \), so that a word is wide enough to hold an index into an input of size \(n\). No instruction is free and none costs more than a constant. Under this model "time" is the number of instructions executed and "space" is the number of cells touched.

Every one of these assumptions is a lie of a useful kind, and knowing which lie is being told is what separates an analysis that predicts behavior from one that does not. Unit-cost arithmetic is false for numbers wider than a word. Multiplying two \(n\)-bit integers is not one operation, which is why Karatsuba's algorithm exists and why the bit complexity of a numerical algorithm can differ from its arithmetic complexity. Unit-cost memory access is false for any machine with a cache hierarchy. On the host used for the measurements on this page a level-one hit is a few cycles and a DRAM miss is on the order of a hundred nanoseconds, a ratio near \(100\times\), which is why an \(O(n \log n)\) algorithm with random access can lose to an \(O(n^2)\) algorithm that streams. The external-memory model of Aggarwal and Vitter (1988) and the cache-oblivious model of Frigo, Leiserson, Prokop, and Ramachandran (1999) repair exactly this by counting block transfers between two levels of memory instead of instructions, and they reproduce the practical superiority of blocked matrix multiply and of merge-based sorting on disk, which the RAM model cannot see. If \(w\) is allowed to be larger than \(\Theta(\log n)\) the model becomes strictly stronger and admits results such as sorting integers in \(O(n \log\log n)\) time (Han, 2002, with Han and Thorup pushing the randomized bound to \(O(n\sqrt{\log\log n})\)), which is not a contradiction of the comparison lower bound because it is not a comparison algorithm. The discipline is to state which model a bound lives in and to check that the model's cost assumptions are the ones the target machine actually charges.

Asymptotic notation, with the quantifiers

Asymptotic notation is a statement about sets of functions, and most confusion about it comes from suppressing the quantifiers. For functions \(f, g : \N \to \R_{\ge 0}\), the definitions are as follows.

Read the difference between \(O\) and \(o\) carefully, because it is the whole point. In \(O\) the constant \(c\) is chosen once by the prover and may be enormous. In \(o\) the bound must hold for every positive \(c\), which is the statement \( \lim_{n\to\infty} f(n)/g(n) = 0 \). So \(3n^2 \in O(n^2)\) and \(3n^2 \notin o(n^2)\), while \(n \log n \in o(n^2)\). The mirrored little-omega \(\omega(g)\) says the ratio diverges. The equality sign in "\(f(n) = O(g(n))\)" is a traditional abuse for \(f \in O(g)\). It is asymmetric and does not compose backwards, so \(O(n) = O(n^2)\) is true read left to right and false read right to left.

A second point trips people. \(O\) is an upper bound, not a claim of tightness, and it is perfectly correct and perfectly useless to say that mergesort runs in \(O(n^{100})\) time. When someone says an algorithm "is \(O(n^2)\) and that is tight" they mean \(\Theta(n^2)\), and when someone says a problem requires \(\Omega(n \log n)\) they are making a much stronger claim, quantified over all algorithms in a model, which is why lower bounds are rare and precious. A third point is that these are statements about the limit, so they say nothing about any particular \(n\). Galactic algorithms make this concrete. The current asymptotically fastest matrix multiplication exponent, around \(\omega < 2.3719\) in the line of work from Coppersmith and Winograd through Alman, Duan, Vassilevska Williams, Xu, Xu, and Zhou (2024), comes with constants so large that no instance humans can build is faster with it than with Strassen or with a well-tuned cubic kernel.

Constants and lower-order terms are dropped because they are properties of a machine and a compiler rather than of an algorithm, so keeping them would make bounds unportable and would not change which algorithm eventually wins. That reasoning fails in three identifiable situations, and each is a real engineering failure mode. When the constant hides a memory-hierarchy effect, two \(O(n \log n)\) sorts can differ by \(5\times\) because one is cache-friendly. When \(n\) is small and stays small, which is why every production sort switches to insertion sort under a threshold near 16 to 32 elements, and why the Karatsuba implementation measured below is slower than schoolbook multiplication until roughly 32 limbs. And when the hidden constant is a function of a parameter treated as constant, an \(O(n)\) algorithm with a \(2^{k}\) factor for a "small" parameter \(k\) is fixed-parameter tractable, and whether that is fast is a question about \(k\), not about \(n\).

Worst case, average case, amortized

Three distinct quantities are all called "the running time" and they answer different questions. Worst case, \( T(n) = \max_{|x| = n} t(x) \), is a guarantee that no input does worse, which is what a service with a latency SLO and an adversarial user population needs. Average case, \( \E_{x \sim \D}[t(x)] \), is an expectation over an assumed input distribution, and it is only as trustworthy as that assumption. The classic cautionary tale is deterministic quicksort with a first-element pivot, average \(\Theta(n \log n)\) under uniformly random permutations and \(\Theta(n^2)\) on the already-sorted inputs that real systems produce constantly.

Amortized analysis is different in kind from both. It is a worst-case bound on a sequence of operations, with no probability anywhere. The claim "push onto a dynamic array is \(O(1)\) amortized" means that any sequence of \(m\) pushes costs \(O(m)\) total, even though individual pushes cost \(\Theta(n)\) when they trigger a resize. The accounting argument runs as follows. On a doubling array, growing from capacity \(2^k\) to \(2^{k+1}\) copies \(2^k\) elements, and the total copy work for \(m\) pushes is \( \sum_{k \le \log_2 m} 2^k < 2m \), so the total is at most \(3m\) unit operations and the average is constant. The three standard techniques are aggregate analysis (bound the total directly, as just done), the accounting method (charge each cheap operation an extra fixed fee and prove the accumulated credit never goes negative, here charging 3 units per push, spending 1 on the write and banking 2, which exactly pays for copying when the array doubles), and the potential method (define \( \Phi \) on the data structure state with \(\Phi \ge \Phi_{\text{initial}}\), and define the amortized cost as \( \hat{c}_i = c_i + \Phi_i - \Phi_{i-1} \), which telescopes to \( \sum \hat{c}_i = \sum c_i + \Phi_n - \Phi_0 \ge \sum c_i \)). For the doubling array take \( \Phi = 2\,(\text{size}) - (\text{capacity}) \). A non-resizing push changes \(\Phi\) by \(+2\) for amortized cost 3, and a resizing push at size \(=\) capacity \(= s\) costs \(s+1\) real units while \(\Phi\) drops from \(s\) to \(2(s+1) - 2s = 2\), giving amortized \( s + 1 + 2 - s = 3 \). Doubling is load-bearing, since growing by a constant additive amount makes the same sequence cost \(\Theta(m^2)\). Fibonacci heaps, splay trees, and union-find with path compression are the classic amortized structures, the last with the inverse-Ackermann bound \( O(\alpha(n)) \) per operation proved by Tarjan (1975), where \(\alpha(n) \le 4\) for every \(n\) that fits in the universe.

Randomized algorithms: Las Vegas and Monte Carlo

A randomized algorithm makes coin flips, and the running time and the output become random variables. The two families are distinguished by which one is allowed to vary. A Las Vegas algorithm is always correct and its running time is random. Randomized quicksort always returns a sorted array and takes \( \Theta(n \log n) \) time in expectation, with a worst case of \(\Theta(n^2)\) that occurs with vanishing probability. A Monte Carlo algorithm has a fixed running time and is allowed to be wrong with bounded probability. Karger's contraction algorithm returns a cut in a fixed \(O(n^2)\) time and that cut is minimum only with probability at least \(2/(n(n-1))\). Monte Carlo algorithms come in one-sided and two-sided flavors. The Miller-Rabin primality test is one-sided (it never calls a prime composite), which is what lets its error be driven to \(2^{-128}\) by independent repetition.

The two families interconvert imperfectly. A Las Vegas algorithm with expected time \(T\) becomes Monte Carlo by truncation. Run for \(kT\) steps and output garbage if it has not finished, which by Markov's inequality fails with probability at most \(1/k\). The reverse conversion needs an efficient verifier. Run the Monte Carlo algorithm, check the answer, and repeat until the check passes, which turns a success probability \(p\) into an expected \(1/p\) repetitions. The important reason to randomize is not that random inputs are easy but that randomization removes the adversary. A deterministic algorithm has a worst-case input that an adversary can construct and a Denial-of-Service attacker can send, while a randomized algorithm's behavior depends on coins the adversary cannot see, so no fixed input is bad. That is the real reason production hash tables use keyed or seeded hash functions, after collision-flooding attacks on unseeded string hashing became a standard exploit against web frameworks around 2011.

Recurrences

Substitution, and the strengthening trick

Divide-and-conquer running times are recurrences, and the most reliable way to solve one is to guess the answer and prove it by induction. The guess is usually obtained from a recursion tree. The induction is where the work is. Consider

Guess \( T(n) \le c\, n \log_2 n \) for all \( n \ge 2 \) and some constant \(c\) to be chosen. The inductive step assumes the bound for all smaller arguments.

where the last inequality holds as soon as \(c \ge 1\). The base case needs care because \(n \log_2 n = 0\) at \(n = 1\), so the bound cannot hold there for any \(c\). The fix is to start the induction at \(n = 2\) and \(n = 3\), which is legitimate because for \(n \ge 4\) the recursion only ever reaches down to those values. With \( T(2) = 2T(1) + 2 = 4 \le c \cdot 2 \cdot 1 \) we need \( c \ge 2 \), and \( T(3) = 2T(1) + 3 = 5 \le c \cdot 3 \log_2 3 = 4.755c \) needs \(c \ge 1.052\). Taking \(c = 2\) satisfies everything, so \( T(n) \le 2 n \log_2 n \).

Now the trick that the method is famous for. Try to prove that \( T(n) = 2T(\lfloor n/2\rfloor) + 1 \) is \(O(n)\) with the guess \(T(n) \le cn\).

The induction fails. The quantity \(cn + 1\) is not \(\le cn\), and no choice of \(c\) repairs it, because the \(+1\) is additive and does not shrink with \(c\). The answer is nonetheless \(\Theta(n)\). The guess is right and the induction hypothesis is too weak. Strengthen it by subtracting a lower-order term, which paradoxically makes the claim harder and the proof easier. Guess \( T(n) \le cn - b \) with \(b > 0\).

whenever \( b \ge 1 \). Now the induction goes through, because the strengthened hypothesis hands back \(2b\) of slack while the recurrence only spends \(b + 1\) of it. The base case \(T(1) = 1 \le c - b\) is satisfied by \(c = 2, b = 1\). This subtract-a-term maneuver is the standard repair whenever an induction on an \(O(\cdot)\) bound fails by an additive constant, and the reason it is needed at all is that big-O notation cannot be carried through an induction. Writing "\(T(n) \le 2 \cdot O(n/2) + 1 = O(n)\)" hides a constant that grows with the depth of the recursion and is a genuine error, not a shortcut.

Recursion trees

A recursion tree is the bookkeeping device that produces the guess. Draw the recurrence as a tree whose root is the non-recursive cost at the top level and whose children are the subproblems. Sum the cost across each level, then sum the levels. For \( T(n) = 3T(n/4) + n^2 \) the tree reads as follows.

level 0:                 n^2                              cost n^2
level 1:      (n/4)^2  (n/4)^2  (n/4)^2                cost 3(n/4)^2 = (3/16) n^2
level 2:   9 nodes of cost (n/16)^2                    cost (3/16)^2 n^2
   ...
level i:   3^i nodes of cost (n/4^i)^2                 cost (3/16)^i n^2
depth:     log_4 n levels, 3^(log_4 n) = n^(log_4 3) leaves

total = n^2 * sum_i (3/16)^i  <  n^2 * 1/(1 - 3/16) = (16/13) n^2  =  Theta(n^2)

The geometric series is dominated by its first term, so the root cost decides and the answer is \(\Theta(n^2)\). Three shapes occur. If the per-level cost decreases geometrically the root dominates and \(T = \Theta(f(n))\). If it increases geometrically the leaves dominate and \(T = \Theta(n^{\log_b a})\), the number of leaves. If it is constant across all \( \log_b n \) levels, every level contributes equally and a factor of \(\log n\) appears. Those three shapes are exactly the three cases of the master theorem, which is nothing more than this argument done once and for all.

The master theorem

Let \( a \ge 1 \), \( b > 1 \) be constants, \(f\) non-negative, and

Write \( n^{\log_b a} \) for the leaf count, sometimes called the watershed function. The three cases follow.

The \(\varepsilon\) in cases 1 and 3 is not decoration. The comparison between \(f\) and \(n^{\log_b a}\) must be polynomial, and when it is merely logarithmic the theorem does not apply. The regularity condition \( a f(n/b) \le c f(n) \) in case 3 says the cost really does shrink going down the tree, which is what licenses summing the geometric series at the root. Applications, all immediate once \(\log_b a\) is computed, fill the table below.

A recurrence the theorem cannot touch is \( T(n) = 2T(n/2) + n/\log n \). Here \(n^{\log_b a} = n\) and \( f(n) = n/\log n \) is smaller than \(n\) but not by any polynomial factor, so case 1 fails. It is not \(\Theta(n \log^k n)\) with \(k \ge 0\), so case 2 fails. A recursion tree settles it. Level \(i\) costs \( 2^i \cdot (n/2^i)/\log(n/2^i) = n/(\log n - i) \), and summing \(i = 0\) to \(\log n - 1\) gives \( n \sum_{j=1}^{\log n} 1/j = \Theta(n \log \log n) \).

Two generalizations are worth naming. Unequal splits, as in \( T(n) = T(n/5) + T(7n/10) + \Theta(n) \) from median-of-medians, are outside the master theorem entirely because the subproblems differ in size. The Akra-Bazzi method (1998) handles \( T(n) = \sum_{i=1}^{k} a_i T(b_i n) + f(n) \) by solving \( \sum_i a_i b_i^{\,p} = 1 \) for the exponent \(p\) and then evaluating \( T(n) = \Theta\!\left( n^p \left( 1 + \int_1^n \frac{f(u)}{u^{p+1}}\,du \right) \right) \). For the median-of-medians recurrence, \( (1/5)^p + (7/10)^p = 1 \) is satisfied at \(p = 1\) exactly when \(1/5 + 7/10 = 9/10 < 1\), so in fact \(p < 1\) is the root and the integral \( \int_1^n u^{-p}\,du = \Theta(n^{1-p}) \) dominates, giving \(T(n) = \Theta(n)\), matching the elementary argument below. Akra-Bazzi also absorbs floors, ceilings, and additive perturbations of the split points, which the master theorem formally does not (the standard remedy there being to prove the floor-free version and then check that floors change nothing, a step usually and reasonably skipped).

Deriving a recurrence from an actual algorithm

Recurrences do not arrive labeled. They are read off a procedure. Take counting inversions in an array, the number of pairs \(i < j\) with \(A[i] > A[j]\), which is the natural measure of how far a ranking is from another and the statistic behind Kendall's tau. The brute force is \(\Theta(n^2)\). The divide-and-conquer version splits the array in half, counts inversions inside the left half and inside the right half recursively, then counts the split inversions with one merge pass. When the merge takes an element from the right half while \(k\) elements remain unconsumed in the left half, those \(k\) elements are each greater than it, so add \(k\) to the count. Every inversion is counted exactly once, because a pair is either inside a half or split across the two, and the merge sees each split pair exactly once. Reading the costs off the procedure gives two recursive calls on inputs of size \(n/2\), plus a merge that touches each element a constant number of times.

by master case 2 with \(\log_b a = \log_2 2 = 1\) and \(k = 0\). The subtle correctness point is that the merge must operate on the sorted halves, so the recursion has to return sorted output as a side effect, which is free because it is mergesort with a counter attached. That coupling of a side effect into a recursive contract is the standard move for turning an \(O(n^2)\) counting problem into an \(O(n \log n)\) one, and the same trick counts points dominated in the plane and computes the number of significant inversions.

Divide and conquer

Karatsuba multiplication

Multiplying two \(n\)-digit integers by the schoolbook method takes \(\Theta(n^2)\) digit products, and in 1956 Kolmogorov conjectured publicly that \(\Omega(n^2)\) was a lower bound. Karatsuba, then a student in Kolmogorov's seminar, refuted the conjecture within about a week. Kolmogorov wrote the result up under both names, and Karatsuba and Ofman (1962) is the citation. The idea is one saved multiplication. Split each operand at the midpoint, \( x = x_1 B + x_0 \) and \( y = y_1 B + y_0 \) where \(B = 10^{n/2}\) or a power of two. The product expands to

which looks like four half-size multiplications, giving \( T(n) = 4T(n/2) + \Theta(n) \), master case 1 with \(\log_2 4 = 2\), so \(\Theta(n^2)\), no progress. The saving is the identity

which recovers the middle coefficient from products already needed plus one new product of half-size numbers. Three recursive multiplications and \(\Theta(n)\) additions give \( T(n) = 3T(n/2) + \Theta(n) \), master case 1 with exponent \( \log_2 3 = 1.58496\ldots \), so \( T(n) = \Theta(n^{1.585}) \). The lesson generalizes. The exponent of a divide-and-conquer algorithm is set by the number of recursive calls, and algebraic identities that trade a multiplication for additions buy polynomial speedups, not constant ones.

The theory was measured on this machine (Intel Xeon Platinum 8480+, single-threaded pure Python, exact integer arithmetic on \(2^{16}\)-base limbs, recursion cutoff 16 limbs). Every product was checked against CPython's own big-integer result.

Fitting \( \log t \) against \( \log n \) over these rows gives a measured exponent of 2.045 for schoolbook against the theoretical 2, and 1.634 for Karatsuba against the theoretical 1.585. The cleanest way to see a derived exponent in measured data is to divide it out. If the theory is right, time divided by the predicted growth function is a constant, so the bars below should be flat even though the raw times span three orders of magnitude.

normalized cost per predicted unit (pure Python, Xeon 8480+); flat = exponent confirmed

  schoolbook,  t / n^2  (ns)              Karatsuba,  t / n^1.585  (ns)
  n=64       82.5  ############           n=64      389.5  ############
  n=128      83.7  ############           n=128     426.5  #############
  n=256      87.3  ############           n=256     436.9  #############
  n=512      89.0  #############          n=512     456.1  ##############
  n=1024     90.4  #############          n=1024    461.2  ##############
  n=2048     97.7  ##############         n=2048    467.5  ##############

  raw time over this range:  schoolbook 0.338 -> 409.6 ms (1212x),  Karatsuba 0.284 -> 82.8 ms (292x)
  normalized cost over the same range:   +18.4%                       +20.0%

The residual drift, the fact that the bars grow by a fifth rather than staying exactly level, is interpreter overhead and allocator pressure that decay slowly relative to \(n\). It is what pushes the fitted exponents to 2.045 and 1.634 rather than the exact 2 and 1.585. The crossover measurements make the small-\(n\) caveat concrete. At 8, 16, and 24 limbs the schoolbook-to-Karatsuba time ratios are 0.971, 0.997, and 0.943, meaning Karatsuba is at best even and sometimes slower. The ratio first exceeds 1 at 32 limbs (1.028) and reaches only 1.283 at 96 limbs. An asymptotically better algorithm with a worse constant needs room to win, which is why GMP switches multiplication algorithms at tuned thresholds rather than using its asymptotically best routine everywhere.

Strassen, and what the measurements show about constants

The same trick applies to matrix multiplication. Splitting \(2n \times 2n\) matrices into four \(n \times n\) blocks, the obvious recursion uses 8 block products, \( T(n) = 8T(n/2) + \Theta(n^2) \), master case 1 with \(\log_2 8 = 3\), the cubic algorithm again. Strassen (1969) found seven bilinear combinations of the blocks whose sums and differences reconstruct all four blocks of the product, giving \( T(n) = 7T(n/2) + \Theta(n^2) \) and \( \Theta(n^{\log_2 7}) = \Theta(n^{2.8074}) \). The identities are unenlightening to stare at. The enlightening fact is that the rank of the \(2\times 2\) matrix multiplication tensor is 7, not 8, and every fast matrix multiplication algorithm since is a statement about tensor rank.

The measurement on this machine separates the two claims hiding in "Strassen is faster." Counting scalar multiplications, the claim is exactly right. At \( n = 256 \) the naive triple loop performs 16,777,216 scalar multiplications and Strassen with a cutoff of 32 performs 11,239,424, a ratio of 1.493, and the fitted exponent of the multiplication count across \(n = 64, 128, 256\) is 2.807, matching \(\log_2 7 = 2.8074\) to three decimals. Counting wall-clock time in pure Python, the claim evaporates, with 955.8 ms naive against 962.3 ms Strassen at \(n = 256\), because the \(\Theta(n^2)\) block additions Strassen adds cost as much per element in an interpreter as the multiplications it saves. Both counts were verified exact against the naive product on integer matrices. In tuned BLAS implementations, where a multiply-add is one fused instruction and additions are pure overhead, Strassen pays off only for \(n\) in the thousands, and the numerically sharper issue is that its subtractions lose relative accuracy on ill-scaled inputs, which is why most production BLAS libraries do not use it by default. The asymptotic record is a separate world. The Coppersmith-Winograd line, continued by Alman and Vassilevska Williams (2021) and by Alman, Duan, Vassilevska Williams, Xu, Xu, and Zhou (2024), currently gives \( \omega < 2.3716 \), with constants that make the algorithms galactic, and DeepMind's AlphaTensor (Fawzi et al., Nature 2022) used reinforcement learning over tensor decompositions to find, among other results, a 47-multiplication algorithm for \(4 \times 4\) matrices over \(\mathbb{F}_2\), one below Strassen's recursive 49.

Selection in worst-case linear time: median of medians

The \(k\)-th smallest of \(n\) elements can be found without sorting. Quickselect partitions around a random pivot and recurses into the one side containing rank \(k\). Its expected time obeys \( T(n) \le T(3n/4) + \Theta(n) \) in the amortized sense that a random pivot lands in the middle half with probability \(1/2\), giving expected \(\Theta(n)\), but an adversarial input drives the deterministic first-element variant to \(\Theta(n^2)\). Blum, Floyd, Pratt, Rivest, and Tarjan (1973) showed the worst case can be made linear by spending linear time choosing a provably good pivot. Split the input into \( \lceil n/5 \rceil \) groups of five, take each group's median by brute force, then recursively select the median of those medians, and use it as the pivot.

The pivot guarantee is combinatorial. At least half of the group medians, about \( \lceil n/5 \rceil / 2 \), are less than or equal to the pivot, and each such group contributes 3 elements at most the pivot (its median and the two below), so at least roughly \( 3n/10 - 6 \) elements are \(\le\) the pivot. Symmetrically at least \(3n/10 - 6\) are \(\ge\) it. The recursion therefore lands on at most \( 7n/10 + 6 \) elements, on top of the recursive call of size \( n/5 \) that found the pivot.

Because \( \tfrac15 + \tfrac{7}{10} = \tfrac{9}{10} < 1 \), the per-level work shrinks geometrically and the total is linear. Substitution makes it precise. Guess \( T(n) \le dn \), then \( T(n) \le d n/5 + 7dn/10 + O(1) + cn = \tfrac{9}{10} d n + cn + O(1) \le dn \) as soon as \( d \ge 10c \) and \(n\) exceeds a constant. Groups of five are the smallest odd size that works. With groups of three the guarantee is only \(n/3\) discarded, the recursion becomes \( T(n/3) + T(2n/3) + cn \) with \( \tfrac13 + \tfrac23 = 1 \), every level of the tree costs \(cn\), and the sum is \(\Theta(n \log n)\).

The measured comparison counts show what the guarantee costs. On random inputs, quickselect used between 7.62 comparisons per element at \( n = 1000 \) and 6.20 at \( n = 64{,}000 \), while median-of-medians used 11.73 rising to 12.95. The deterministic algorithm pays roughly a factor of two on inputs where randomization was already safe. That ratio is the entire story of why practical libraries use introselect, quickselect with a median-of-medians fallback triggered by slow progress, which keeps the expected constant of the randomized method and the worst-case linearity of BFPRT. C++ std::nth_element is specified to be exactly this shape.

The comparison lower bound, and a measurement that touches it

Sorting has a genuine \(\Omega(n \log n)\) lower bound in the comparison model, and the proof is worth internalizing because true lower bounds are rare. Any deterministic comparison sort, run on inputs that are permutations of \( \{1,\dots,n\} \), can be drawn as a binary decision tree, with internal nodes the comparisons and leaves the sorted orders the algorithm commits to. Two different permutations must reach different leaves, since otherwise the algorithm performs the same rearrangement on both and gets one of them wrong. So the tree has at least \( n! \) leaves, and a binary tree of height \(h\) has at most \(2^h\) leaves, forcing

the middle equality by Stirling. The worst-case number of comparisons is the tree height, so it is at least \( \log_2 n! \), and an averaging argument over leaves gives the same bound for the average case. The bound is about the model. It counts only two-way branches on comparisons, so radix sort and Han's \(O(n \log\log n)\) integer sorting do not contradict it, they refuse the model.

The measured comparison counts sit close to this floor. At \( n = 16{,}000 \), \( \log_2 n! = 200{,}377.7 \), and top-down mergesort averaged 203,309.1 comparisons over 25 random permutations, a ratio of 1.0146. The smooth upper bound \( n\log_2 n - n + 1 \) evaluates to 207,453.5 there. Randomized quicksort averaged 262,497.1 comparisons, and the exact expectation derived below, \( 2(n+1)H_n - 4n = 264{,}263.4 \) at this \(n\), is within 0.7% of the measurement. Quicksort runs at about 1.31 times the decision-tree floor at these sizes, drifting toward its asymptotic \( 2\ln 2 \approx 1.386 \). Wall clock tells the constants-versus-model story one more time. Normalizing the measured comparison counts by the information-theoretic floor \( \log_2 n! \) makes the two algorithms' distance from optimality visible as a flat line, which is the empirical content of the lower bound.

measured comparisons divided by the decision-tree floor log2(n!), 25 random permutations each

  n=1000    mergesort 1.0214 ####################   quicksort 1.3131 ##########################
  n=2000    mergesort 1.0189 ####################   quicksort 1.3343 ###########################
  n=4000    mergesort 1.0173 ####################   quicksort 1.3090 ##########################
  n=8000    mergesort 1.0157 ####################   quicksort 1.3164 ##########################
  n=16000   mergesort 1.0146 ####################   quicksort 1.3100 ##########################

  mergesort creeps toward 1.0; quicksort sits at its asymptotic 2 ln 2 = 1.386 constant, less
  the O(n) correction in 2(n+1)H_n - 4n, which is why the ratio is near 1.31 at these sizes

No algorithm in the model can push the left column below 1.0, and mergesort's steady approach to it is why the \( \Omega(n \log n) \) bound is called tight. Wall clock says something orthogonal. At \( n = 160{,}000 \) the pure-Python mergesort took 546.18 ms, quicksort 441.72 ms, and the built-in Timsort 27.57 ms, a \(20\times\) gap at identical asymptotics, from doing the same comparisons in C and exploiting pre-sorted runs.

The fast Fourier transform

The discrete Fourier transform of a coefficient vector \( a = (a_0, \dots, a_{n-1}) \) is its evaluation at the \(n\)-th roots of unity \( \omega_n^k = e^{2\pi i k / n} \).

Evaluated naively this is a dense \(n \times n\) matrix-vector product, \(\Theta(n^2)\). The Cooley-Tukey observation (1965, and in Gauss's unpublished notes of 1805) is that the roots of unity have exactly the recursive structure a divide-and-conquer algorithm needs. Split the polynomial \( A(x) = \sum_j a_j x^j \) by parity of index into \( A_e(y) = a_0 + a_2 y + a_4 y^2 + \cdots \) and \( A_o(y) = a_1 + a_3 y + \cdots \), so that

The halving lemma does the rest. The squares of the \(n\)-th roots of unity are precisely the \(n/2\)-th roots of unity, each hit twice, because \( (\omega_n^{k + n/2})^2 = \omega_n^{2k} \omega_n^{n} = \omega_n^{2k} \). So evaluating \(A\) at all \(n\) roots needs \(A_e\) and \(A_o\) evaluated at only \(n/2\) points, two half-size subproblems, plus \(n\) multiply-adds to combine, the butterfly.

using \( \omega_n^{k+n/2} = -\omega_n^k \). The recurrence is \( T(n) = 2T(n/2) + \Theta(n) \), master case 2, \( \Theta(n \log n) \). The inverse transform is the same algorithm with conjugated roots and a final division by \(n\), which follows from the orthogonality relation below.

where the second case is the geometric series formula and \( \omega_n^{n(k-k')} = 1 \). In matrix language, the DFT matrix \(F\) satisfies \( F^* F = nI \), so \( F^{-1} = F^*/n \). The payoff is the convolution theorem. Multiplying two polynomials is pointwise multiplication of their evaluations, so polynomial (and hence big-integer) multiplication costs two forward FFTs, \(n\) pointwise products, and one inverse FFT, \( \Theta(n \log n) \) against the schoolbook \(\Theta(n^2)\) and Karatsuba's \(\Theta(n^{1.585})\).

Floating-point FFT multiplication has a correctness obligation the other algorithms do not. Coefficients come back as doubles and are rounded to integers, which is valid only when accumulated rounding error stays below one half. The measured run used 12-bit limbs precisely so that every pointwise product stays under the \(2^{53}\) mantissa limit. Across 1,536 to 24,576-bit products, the worst observed rounding error was \(4.65 \times 10^{-4}\), and every rounded result matched CPython's exact integers, as did a 49,152-bit stress test with maximum error \(2.63 \times 10^{-4}\). The fitted time exponent was 1.138 against the theoretical \( n \log n \) (exponent 1 plus a logarithmic factor). The head-to-head at equal bit widths is the cleanest summary of this whole section. Multiplying two 24,576-bit integers took 372.16 ms schoolbook, 79.69 ms Karatsuba, 17.78 ms FFT in pure Python, and 2.92 ms with CPython's built-in big integers, which run Karatsuba in C. Two different asymptotic classes and one constant-factor class are all visible in a single row. The asymptotic endpoint of integer multiplication is Harvey and van der Hoeven's \( O(n \log n) \) algorithm (Annals of Mathematics, 2021), galactic for now, closing a conjecture of Schönhage and Strassen from 1971 whose own \( O(n \log n \log\log n) \) algorithm is what GMP actually runs above its FFT threshold.

Randomization as a design resource

Randomized quicksort, analyzed with indicator variables

The expected comparison count of quicksort with uniformly random pivots has an exact closed form, and the derivation is the canonical example of the indicator-variable method. Let \( z_1 < z_2 < \cdots < z_n \) be the input elements in sorted order and let \( X_{ij} \) indicate that \( z_i \) and \( z_j \) are ever compared. Two elements are compared exactly when one of them is chosen as a pivot while both are in the same subarray, and \(z_i, z_j\) stay in the same subarray precisely until the first pivot chosen from the block \( Z_{ij} = \{ z_i, z_{i+1}, \dots, z_j \} \). That first pivot is uniform over the \( j - i + 1 \) elements of \(Z_{ij}\). The pair is compared if it is \(z_i\) or \(z_j\) and separated without comparison otherwise. Hence

where the last equality is bookkeeping with the harmonic number \( H_n = \sum_{k \le n} 1/k \), and since \( H_n = \ln n + \gamma + O(1/n) \), the expectation is \( 2n \ln n - O(n) \approx 1.386\, n \log_2 n \). The measurement agrees closely. Averaged over 25 random permutations, the measured counts at \( n = 1000, 4000, 16000 \) were 11,199.7, 55,107.8, and 262,497.1 against formula values 10,985.9, 54,988.9, and 264,263.4, ratios 1.020, 1.002, and 0.993. An asymptotic analysis with the constant kept is a falsifiable prediction, and this one survives contact with the machine.

Universal hashing, with the collision bound proved

Hash tables are the other place randomization buys a worst-case guarantee, and the guarantee is about the family of hash functions rather than any one of them. No fixed \( h : U \to \{0,\dots,m-1\} \) is safe when \( |U| > m \). The pigeonhole principle puts at least \( |U|/m \) keys in some bucket, and an adversary who knows \(h\) sends exactly those. A family \( \mathcal{H} \) is universal (Carter and Wegman, 1979) if for every pair of distinct keys \( x \ne y \),

which is what a uniformly random function would give. The canonical family is \( h_{a,b}(x) = ((ax + b) \bmod p) \bmod m \) with \(p\) a prime exceeding \( \max U \), \( a \in \{1,\dots,p-1\} \), \( b \in \{0,\dots,p-1\} \), and it is universal for an exact reason. Fix \( x \ne y \) and write \( r = (ax+b) \bmod p \), \( s = (ay+b) \bmod p \). Then \( r - s \equiv a(x - y) \pmod p \), which is non-zero because \( \mathbb{Z}_p \) is a field, \( a \ne 0 \), and \( x \ne y \), so \( r \ne s \) always and no collision happens before the final mod. The map \( (a,b) \mapsto (r,s) \) is a bijection from the \( p(p-1) \) parameter choices onto the \( p(p-1) \) ordered pairs of distinct residues (given \(r, s\) solve back for \(a = (r-s)(x-y)^{-1}\) and \(b = r - ax\)), so \( (r,s) \) is uniform over distinct pairs. A collision requires \( r \equiv s \pmod m \). For each of the \(p\) values of \(r\), the number of \( s \ne r \) with \( s \equiv r \pmod m \) is at most \( \lceil p/m \rceil - 1 \le (p-1)/m \). Dividing by the \( p - 1 \) equally likely values of \(s\) gives probability at most \( 1/m \).

The consequence is the bound that makes chained hashing work. For \(n\) keys, let \( X_{xy} \) indicate that \(x\) and \(y\) collide. Then \( \E[\#\text{colliding pairs}] = \sum_{x < y} \P(h(x) = h(y)) \le \binom{n}{2}/m \), and the expected number of keys sharing a bucket with a given key \(x\) is at most \( (n-1)/m < \alpha \), the load factor, so an expected search costs \( \Theta(1 + \alpha) \) with no assumption about the input distribution. The measurement is a direct test of the inequality. Placing 500 keys into 64 buckets over 2,000 independent draws of \( (a,b) \) with \( p = 2^{61} - 1 \) gave a mean of 1,948.6 colliding pairs against the bound \( \binom{500}{2}/64 = 1{,}949.2 \), a ratio of 0.9997. The family is not merely universal but essentially exactly so, which is why it, or a faster relative such as SipHash or a multiply-shift scheme, sits under every hash table that has to survive untrusted keys.

Karger's contraction algorithm

Karger's algorithm (1993) finds a global minimum cut of an undirected multigraph by doing something that sounds too crude to work. Repeatedly pick a uniformly random edge and contract its endpoints into one supernode (keeping parallel edges, discarding self-loops) until two supernodes remain, then output the edges between them, which form a cut of the original graph. Fix any particular minimum cut \(C\) with \( |C| = k \). The algorithm outputs \(C\) exactly when no edge of \(C\) is ever contracted. When \(m_i\) edges and \( n - i \) supernodes remain, every supernode still has degree at least \(k\) (contracting cannot reduce the min cut below \(k\)), so \( m_i \ge k(n-i)/2 \), and the chance of hitting \(C\) at step \(i\) is at most \( k / m_i \le 2/(n-i) \). Surviving all \( n - 2 \) contractions therefore has probability at least

a telescoping product in which every numerator cancels a denominator two steps later. Two corollaries fall out. Repetition amplifies, since \( t = \binom{n}{2} \ln(1/\delta) \) independent runs miss a fixed min cut with probability at most \( (1 - 2/(n(n-1)))^t \le e^{-\ln(1/\delta)} = \delta \). And counting comes free. Each distinct minimum cut is output with probability at least \( 2/(n(n-1)) \), and these events are disjoint, so a graph has at most \( \binom{n}{2} \) minimum cuts, a purely combinatorial fact proved by a randomized algorithm's analysis.

The measured behavior on a 10-vertex, 24-edge graph with min cut 2 bears this out. Over 20,000 independent runs the empirical per-run success probability was 0.57185, more than 25 times the proved floor of \( 2/(10 \cdot 9) = 0.02222 \), which is the usual relationship between a worst-case bound and a typical instance. Running the prescribed 208 repetitions for a 1% failure target predicts failure probability \( (1 - 0.02222)^{208} = 0.00933 \), and the observed failure rate was 0. The cycle \(C_{10}\) shows the bound is not slack in general. It has \( \binom{10}{2} = 45 \) distinct minimum cuts (any two arc-separating vertex pairs), and over 50,000 runs the algorithm returned some min cut every single time, but any one designated cut only 2.21% of the time, essentially exactly the \( 2/(n(n-1)) = 2.22\% \) floor. The bound is tight per cut, and the cycle saturates the \( \binom{n}{2} \) count. The fastest descendant of this idea, Karger-Stein recursive contraction, reuses the observation that early contractions are safe and late ones are risky, re-branching only near the end for \( O(n^2 \log^3 n) \) total time.

Greedy algorithms and when they can be trusted

Interval scheduling and the exchange argument

A greedy algorithm commits to a locally best choice and never reconsiders. Most greedy ideas are wrong, so the subject is really the proof techniques that certify the rare correct ones. The cleanest specimen is interval scheduling. Given \(n\) intervals \( [s_i, f_i) \), select a maximum-size subset that is pairwise disjoint. The correct greedy rule is to sort by finish time and repeatedly take the earliest-finishing interval compatible with what has been taken. Three plausible rules fail, and knowing the counterexamples is as important as knowing the proof. Earliest start time fails against one long interval that starts first and blocks everything. Shortest interval fails against a short interval straddling the boundary of two long compatible ones. Fewest conflicts fails against a four-layer construction in which the least-conflicted interval sits in the middle of the unique optimal chain.

The correctness proof is a greedy stays ahead induction, one of the two standard templates. Let the greedy choices in order be \( g_1, \dots, g_k \) with finish times \( f(g_1) \le \cdots \le f(g_k) \), and let \( o_1, \dots, o_m \) be any optimal solution sorted by finish time. The claim, proved by induction on \(r\), is that \( f(g_r) \le f(o_r) \) for every \( r \le \min(k, m) \). For \( r = 1 \), greedy takes the globally earliest finisher. For the step, assume \( f(g_{r-1}) \le f(o_{r-1}) \). Then \( o_r \) starts at or after \( f(o_{r-1}) \ge f(g_{r-1}) \), so \( o_r \) was compatible with greedy's first \( r-1 \) picks at the moment greedy made its \(r\)-th choice. Greedy chose the earliest-finishing compatible interval, so \( f(g_r) \le f(o_r) \). Now suppose \( m > k \). The interval \( o_{k+1} \) starts at or after \( f(o_k) \ge f(g_k) \), so it was compatible with everything greedy took, and greedy, which only halts when nothing compatible remains, would have taken it, a contradiction. Hence \( k = m \) and greedy is optimal. The proof's engine is an invariant comparing greedy's prefix to an arbitrary optimum's prefix, resource by resource. The same skeleton proves the deadline-scheduling and fractional-knapsack greedies.

The second template is the exchange argument. Take any optimal solution that differs from greedy's, find the first difference, and swap the optimum's choice for greedy's without making it worse, concluding after finitely many swaps that some optimum agrees with greedy everywhere. For minimizing maximum lateness by earliest-deadline-first (jobs with processing times and deadlines on one machine), the exchange is between an adjacent inverted pair. If the schedule runs job \(j\) immediately before job \(i\) with \( d_i < d_j \), swapping them cannot increase any completion time except \(j\)'s new completion, which equals \(i\)'s old completion, and \( C_{\text{old}}(i) - d_i \ge C_{\text{new}}(j) - d_j \) since \( d_j > d_i \) and \( C_{\text{new}}(j) = C_{\text{old}}(i) \). So the maximum lateness does not increase. Bubble-sorting the optimum into deadline order with such swaps proves EDF optimal. Exchange arguments are the more general tool, and they are exactly what matroid theory axiomatizes.

Huffman coding

Huffman's algorithm (1952) builds a minimum-expected-length prefix code by repeatedly merging the two least frequent symbols into a combined pseudo-symbol. Its correctness proof is the exchange argument in two lemmas. Lemma 1 (greedy choice). Some optimal tree has the two rarest symbols \(x, y\) as sibling leaves at maximum depth. Take any optimal tree and let \(a, b\) be siblings at maximum depth. Swapping \(x\) with \(a\) changes the cost by \( (p_a - p_x)(d_x - d_a) \le 0 \), since \( p_x \le p_a \) and \( d_a \ge d_x \), and likewise for \(y\) with \(b\). So the swaps produce an optimal tree of the desired form. Lemma 2 (optimal substructure). If \(T\) is optimal for the alphabet with \(x, y\) merged into \(z\) with \( p_z = p_x + p_y \), then expanding \(z\) back into siblings \(x, y\) gives an optimal tree for the original alphabet, because the two costs differ by exactly the constant \( p_x + p_y \), as \( B(T) = B(T') + p_x + p_y \), so minimizing one minimizes the other. Induction on alphabet size completes the proof. The information-theoretic frame is Shannon's source coding theorem, which lower-bounds any prefix code's expected length by the entropy \( H(p) = -\sum_i p_i \log_2 p_i \), and Huffman achieves within 1 bit of it, \( H(p) \le L_{\text{Huffman}} < H(p) + 1 \), the slack coming from codeword lengths being forced to integers. Arithmetic coding removes that constraint and gets within any \(\varepsilon\).

Matroids: the exact boundary of greedy correctness

There is a theorem that says precisely when the generic greedy algorithm works, and it is one of the deeper facts in the subject. A matroid is a pair \( M = (E, \mathcal{I}) \) of a finite ground set and a family of independent subsets satisfying three axioms. First, \( \emptyset \in \mathcal{I} \). Second, heredity, meaning subsets of independent sets are independent. Third, exchange, meaning that if \( A, B \in \mathcal{I} \) and \( |A| < |B| \), some \( x \in B \setminus A \) has \( A \cup \{x\} \in \mathcal{I} \). The two motivating examples are linearly independent column subsets of a matrix (exchange is Steinitz exchange from linear algebra) and acyclic edge subsets of a graph (exchange holds because a forest with fewer edges has fewer components, so some edge of the larger forest bridges two components of the smaller). The exchange axiom forces all maximal independent sets, the bases, to have equal size, mirroring dimension.

Theorem (Rado 1957, Edmonds 1971). For a hereditary system \( (E, \mathcal{I}) \), the greedy algorithm (sort by weight descending, add each element if independence is preserved) returns a maximum-weight independent set for every non-negative weighting if and only if \( (E, \mathcal{I}) \) is a matroid.

If. Let greedy pick \( g_1, \dots, g_k \) in weight order \( w(g_1) \ge \cdots \ge w(g_k) \) and let \( B = \{o_1, \dots, o_k\} \) be any max-weight basis, also sorted descending. Claim \( w(g_r) \ge w(o_r) \) for all \(r\). If not, take the least \(r\) with \( w(g_r) < w(o_r) \), and apply exchange to \( A = \{g_1, \dots, g_{r-1}\} \) and \( B' = \{o_1, \dots, o_r\} \). Some \( o_j \in B' \setminus A \) keeps \( A \cup \{o_j\} \) independent, and \( w(o_j) \ge w(o_r) > w(g_r) \), so greedy, scanning in weight order, would have accepted \( o_j \) before ever reaching \( g_r \), a contradiction. Summing the termwise inequalities, greedy's weight is at least the optimum's. Only if. Suppose exchange fails for some \( A, B \in \mathcal{I} \) with \( |A| = p < |B| \) and no element of \( B \setminus A \) extendable into \(A\). Rig the weights, \( w = 1 + \tfrac{1}{2p} \) on \(A\), \( w = 1 \) on \( B \setminus A \), and \( w = 0 \) elsewhere. Greedy takes all of \(A\) first, then can add nothing of \( B \setminus A \), finishing with at most \( p (1 + \tfrac1{2p}) = p + \tfrac12 \). But \(B\) alone weighs at least \( |B| \ge p + 1 \). Greedy loses, so correctness for all weights forces the exchange axiom.

The theorem explains the landscape. Spanning trees form a matroid (the graphic matroid), so Kruskal's algorithm is an instance of matroid greedy and inherits correctness. Interval scheduling does not form a matroid (two disjoint short intervals and one long conflicting interval violate exchange), which is why its greedy needs the bespoke stays-ahead proof and why the greedy rule had to be chosen so carefully. Matchings in bipartite graphs are not a matroid either, but they are the intersection of two matroids, and Edmonds proved max-weight matroid intersection is still polynomial, a two-matroid greedy. Three-matroid intersection is NP-hard, and that is exactly where TSP lives (Hamiltonian path is a three-matroid intersection). The submodular story continues the thread. Greedy on a monotone submodular function over a matroid constraint gives a \( (1 - 1/e) \)-approximation (Nemhauser, Wolsey, Fisher, 1978), the result behind sensor placement and influence maximization, and Fujishige and Schrijver's polymatroid theory generalizes all of it.

Dynamic programming

Optimal substructure, derived rather than assumed

Dynamic programming applies when a problem's optimal solution is built from optimal solutions to a polynomial-size family of subproblems. The discipline is a four-step derivation, and every DP below follows it. First, define the subproblem in words, with explicit parameters. Second, write the recurrence by conditioning on the last decision, arguing that each branch reduces to a strictly smaller subproblem of the same family (this cut-and-paste argument is the optimal substructure proof, since if the remainder were not optimal, splicing in a better remainder would improve the whole, a contradiction). Third, identify base cases. Fourth, order the computation so dependencies are ready, either bottom-up or by memoized recursion, and count states times transition cost. When the second step fails, DP fails. Longest simple path has no such decomposition, because the subpaths of an optimal simple path constrain each other through the vertices they consume, and indeed longest path is NP-hard while shortest path, whose subpaths are free of such interaction (when no negative cycles exist), is the canonical DP.

Fibonacci shows why memoization matters at all. The naive recursion recomputes \( F(n-2) \) from both branches and costs \( \Theta(\varphi^n) \), while the same recurrence over \(n\) distinct states costs \( \Theta(n) \). The exponential-to- polynomial collapse is always the same phenomenon, a recursion tree with exponentially many nodes but only polynomially many distinct subproblems.

Weighted interval scheduling

Add weights to interval scheduling and greedy dies immediately (a heavy interval finishing late beats two light early ones). Sort intervals by finish time and let \( p(j) \) be the largest index \( i < j \) with \( f_i \le s_j \), the latest interval compatible with \(j\). Define \( \text{OPT}(j) \) as the best achievable value using only intervals \( 1..j \). Condition on the last decision, whether interval \(j\) is used.

If \(j\) is taken, nothing between \( p(j)+1 \) and \( j-1 \) can be, and the rest of the solution is an optimal solution over \( 1..p(j) \) by cut-and-paste. With \( p(\cdot) \) computed by binary search after sorting, the total is \( O(n \log n) \). Reconstructing the actual set is a backward walk. At \(j\), compare the two branches and descend into whichever achieved the max. This two-line recurrence is the template for a large family. The "take it or leave it, and taking it jumps you backward" shape reappears in knapsack, in RNA folding, and in every segmentation DP.

Longest common subsequence, and reading the answer back out

A subsequence deletes characters without reordering the rest. For strings \(a\) of length \(m\) and \(b\) of length \(n\), let \( L[i][j] \) be the length of the longest common subsequence of the prefixes \( a_{1..i} \) and \( b_{1..j} \). Condition on the last characters. If \( a_i = b_j \), some optimal LCS uses that matched pair, for if it did not, appending the pair to it would give a longer common subsequence of the same prefixes, a contradiction, so \( L[i][j] = L[i-1][j-1] + 1 \). If \( a_i \ne b_j \), the optimal subsequence cannot use both characters as its last symbol, so at least one of them is unused and the answer is the better of dropping either.

\( O(mn) \) time and space. The traceback is a separate algorithm run on the finished table, walking from \( (m,n) \) toward the origin. At a cell with \( a_i = b_j \), emit \(a_i\) and step diagonally. Otherwise step to whichever of \( L[i-1][j] \), \( L[i][j-1] \) achieved the max. The walk takes \( O(m+n) \) steps because each step decreases \( i + j \). Note the asymmetry with the DP itself. The table can be computed in \( O(\min(m,n)) \) space by keeping one row, but then the traceback is impossible, which is exactly the trade Hirschberg's algorithm resolves by recursively locating the crossing point of the optimal path through the middle column, paying a factor of two in time for linear space. For \( a = \) AGGTAB and \( b = \) GXTXAYB, the table's final entry is 4 and the traceback emits GTAB. The implementation was checked against brute-force enumeration of all \( 2^{m} \) subsequences on 100 random small instances with zero mismatches. LCS on \(k\) sequences is \( O(n^k) \) and NP-hard when \(k\) is part of the input, and the length of the LCS of two permutations is where patience sorting gets \( O(n \log n) \) by turning the problem into longest increasing subsequence.

Edit distance

The Levenshtein distance between strings \(a\) (length \(m\)) and \(b\) (length \(n\)) is the minimum number of single-character insertions, deletions, and substitutions transforming one into the other. The subproblem is \( D[i][j] \), the distance between the prefixes \( a_{1..i} \) and \( b_{1..j} \). Condition on how an optimal edit script treats the last characters. It either deletes \( a_i \) (cost 1, reducing to \( D[i-1][j] \)), inserts \( b_j \) (cost 1, reducing to \( D[i][j-1] \)), or aligns \( a_i \) with \( b_j \) (cost 0 if equal, else 1, reducing to \( D[i-1][j-1] \)).

The state space is \( (m+1)(n+1) \) and each transition is \(O(1)\), so \( O(mn) \) time. Keeping only the previous row gives \( O(\min(m,n)) \) space, and Hirschberg's divide-and-conquer trick recovers the full alignment in linear space by finding where the optimal path crosses the middle row. The computation for kitten \(\to\) sitting run on this machine produces the table below (rows are prefixes of kitten, columns of sitting), with final answer \( D[6][7] = 3 \), realized by substituting k\(\to\)s and e\(\to\)i and inserting g. The same code gives distance 5 for intention \(\to\) execution, the classic alignment example, and the implementation was verified against brute-force edit enumeration on 100 random instances with zero mismatches.

        ""  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   <- D[6][7] = 3

Each cell visibly obeys the recurrence. For instance \( D[5][5] \) (prefixes kitte, sitti) is \( \min(D[4][5]{+}1, D[5][4]{+}1, D[4][4]{+}1) = \min(3, 3, 2) = 2 \), the diagonal branch paying 1 because e \(\ne\) i. Whether \(O(mn)\) can be beaten is a story told in the fine-grained section below. Backurs and Indyk proved that a strongly subquadratic edit distance algorithm would refute the Strong Exponential Time Hypothesis.

0/1 knapsack and pseudo-polynomial time

Given \(n\) items with weights \( w_i \) and values \( v_i \) and capacity \(W\), maximize value subject to total weight at most \(W\), each item used at most once. The subproblem \( K[i][c] \) is the best value using items \( 1..i \) within capacity \(c\). Condition on item \(i\).

Time is \( O(nW) \), which looks polynomial and is not. The input encodes \(W\) in \( \log W \) bits, so \( O(nW) \) is exponential in the input size, a pseudo-polynomial bound. This is not pedantry. Knapsack is NP-hard, and the DP is consistent with that precisely because of the encoding distinction. Problems like knapsack that are hard only when numbers are huge are weakly NP-hard, and they admit an FPTAS. Scaling values by \( \mu = \varepsilon v_{\max} / n \), rounding down, and running the value-indexed DP \( O(n^2 v_{\max}/\mu) = O(n^3/\varepsilon) \) loses at most \( n\mu = \varepsilon v_{\max} \le \varepsilon \cdot \text{OPT} \), giving a \( (1-\varepsilon) \)-approximation in time polynomial in \(n\) and \( 1/\varepsilon \). Strongly NP-hard problems (TSP, bin packing in full generality) provably have no FPTAS unless P \(=\) NP, because an FPTAS with \( \varepsilon \) below the resolution of the objective would solve them exactly in polynomial time. The knapsack DP here was verified against exhaustive subset enumeration on 200 random instances with zero mismatches.

Matrix chain order: a DP over intervals

Matrix multiplication is associative, so the product \( A_1 A_2 \cdots A_n \) with \( A_i \) of shape \( p_{i-1} \times p_i \) can be parenthesized in \( C_{n-1} \) (Catalan-many) ways, all giving the same matrix at very different costs, since multiplying a \( p \times q \) by a \( q \times r \) matrix costs \( pqr \) scalar multiplications. The arithmetic is stark. For shapes \( 10 \times 100 \), \( 100 \times 5 \), \( 5 \times 50 \), the order \( ((A_1 A_2) A_3) \) costs \( 10 \cdot 100 \cdot 5 + 10 \cdot 5 \cdot 50 = 5000 + 2500 = 7{,}500 \) while \( (A_1 (A_2 A_3)) \) costs \( 100 \cdot 5 \cdot 50 + 10 \cdot 100 \cdot 50 = 25{,}000 + 50{,}000 = 75{,}000 \), a factor of ten from parentheses alone. Condition on the outermost multiplication, the last one performed, which splits the chain at some \(k\).

Optimal substructure holds by cut-and-paste. The two sides of the split are independent subchains, and replacing either with a cheaper parenthesization would improve the whole. The state is an interval, so this is the interval-DP template rather than the prefix template used above. The computation order must be by increasing chain length \( j - i \), since \( m[i][j] \) depends on strictly shorter intervals. There are \( \Theta(n^2) \) states and \( O(n) \) transitions each, so \( \Theta(n^3) \) time and \( \Theta(n^2) \) space, against \( \Theta(4^n / n^{1.5}) \) parenthesizations. The same interval shape solves optimal binary search trees, RNA secondary structure folding, and the CYK parsing algorithm for context-free grammars. Recognizing "the last operation splits the range in two" is the whole content of spotting it.

Held-Karp: exponential DP done honestly

TSP on \(n\) cities has \( (n-1)! \) tours. The Bellman-Held-Karp dynamic program (1962) cuts this to \( O(2^n n^2) \), still exponential but a different exponential. At \( n = 20 \), \( 19! \approx 1.2 \times 10^{17} \) against \( 2^{20} \cdot 400 \approx 4.2 \times 10^8 \). Fix city 0 as the start. The subproblem \( g(S, j) \), for \( 0 \in S \subseteq \{0..n{-}1\} \) and \( j \in S \), is the cheapest path that starts at 0, visits exactly the cities of \(S\), and ends at \(j\). Condition on the city visited just before \(j\).

The cut-and-paste argument holds because the path's prefix before \(j\) must itself be a cheapest path over \( S \setminus \{j\} \) ending at \(k\). Crucially the state records only which set was visited and the endpoint, not the order, and that forgetting is where the factorial collapses to \(2^n\). There are \( 2^n n \) states and \( O(n) \) transitions each. The measured scaling on this machine confirms the shape. Dividing wall-clock time by \( 2^n n^2 \) units gives 30.5, 29.2, 28.7, 28.5, 28.4, 29.0 nanoseconds-scale units per state-transition for \( n = 10, 12, 13, 14, 15, 16 \), flat as theory demands, with absolute times growing from 3.13 ms at \( n = 10 \) to 486.83 ms at \( n = 16 \). Each added city roughly doubles the time, which is what \( 2^n \) means operationally. The bitmask DP was verified against brute-force permutation search on 200 random instances with zero mismatches. No algorithm for exact TSP is known that beats \( O(2^n \text{poly}(n)) \) in general, and whether \( 1.9999^n \) is achievable is open.

Graph search and the structure it exposes

BFS, DFS, and the edge classification

Breadth-first search from \(s\) processes vertices in a FIFO queue and assigns \( d[v] = d[u] + 1 \) when it first reaches \(v\) from \(u\). Claim. \( d[v] \) equals the hop distance \( \delta(s,v) \). One direction is immediate, since \( d[v] \) is realized by an actual path. For the other, induct on \( \delta(s,v) \). A shortest path to \(v\) passes through some \(u\) with \( \delta(s,u) = \delta(s,v) - 1 \), which by induction was enqueued with \( d[u] = \delta(s,u) \). The queue is monotone (keys are dequeued in non-decreasing \(d\), because each insertion appends a key one larger than the key being processed), so \(v\) is discovered no later than when \(u\) is scanned, hence \( d[v] \le d[u] + 1 = \delta(s,v) \). That monotone-queue argument is the unweighted special case of Dijkstra's invariant, and it degrades gracefully. With weights in \( \{0,1\} \) a deque replaces the queue (0-1 BFS), and with small integer weights a bucket queue does, which is Dial's algorithm.

Depth-first search timestamps each vertex with a discovery time \( d[v] \) and a finish time \( f[v] \). The parenthesis theorem says that for any \(u,v\) the intervals \( [d[u], f[u]] \) and \( [d[v], f[v]] \) are either disjoint or nested, never crossing, because the recursion is a stack discipline. Nesting means descendancy in the DFS forest, and that single fact classifies every edge \( (u,v) \) by the color of \(v\) when the edge is explored. White gives a tree edge, gray a back edge (to an ancestor, since the gray vertices are exactly the current stack), and black either a forward edge (to a descendant, \( d[u] < d[v] \)) or a cross edge (\( d[v] < d[u] \), a different subtree already finished). In an undirected graph only tree and back edges occur, since either endpoint would have explored the edge first from whichever side it reached. The immediate corollary is the standard cycle test, that a directed graph has a cycle if and only if a DFS produces a back edge. Necessity is immediate. Sufficiency follows because the first vertex of the cycle to be discovered stays gray until all others on the cycle finish, so the cycle's edge into it is a back edge.

Topological order and strongly connected components

In a DAG, listing vertices in decreasing finish time is a topological order. For the proof, take any edge \( (u,v) \) and consider its color case. \(v\) white means \(v\) becomes a descendant and finishes inside \(u\)'s interval, so \( f[v] < f[u] \). \(v\) gray is impossible, since it would be a back edge and the graph is acyclic. \(v\) black means \(v\) already finished, so again \( f[v] < f[u] \). Every edge therefore points from later finish to earlier, which is exactly the definition. Kahn's algorithm, repeatedly emitting a vertex of in-degree zero, is the same theorem in queue form and detects cycles by ending early.

Strongly connected components are the equivalence classes of mutual reachability, and the condensation, one vertex per component, is always a DAG (a cycle among components would merge them). Kosaraju's algorithm runs DFS on \(G\) recording finish times, then runs DFS on the reverse graph \( G^{\mathsf{T}} \) taking roots in decreasing finish order. Each resulting tree is one SCC. Correctness rests on one lemma. If \(C\) and \(C'\) are distinct SCCs with an edge from \(C\) to \(C'\), then \( \max_{v \in C} f[v] > \max_{v \in C'} f[v] \). In case one, the first vertex discovered among \( C \cup C' \) lies in \(C\). Then every vertex of \(C'\) is reachable from it while it is gray, so all of \( C \cup C' \) finishes within its interval and the maximum is in \(C\). In case two, the first lies in \(C'\). Since the condensation is acyclic there is no path back from \(C'\) to \(C\), so that DFS finishes all of \(C'\) without touching \(C\), and \(C\) finishes strictly later. Given the lemma, the highest-finishing vertex overall lies in a source component of the condensation, which is a sink in \( G^{\mathsf{T}} \), so the second DFS from it reaches exactly that component and stops. Peeling components in that order and induction finish the proof. Tarjan's algorithm (1972) gets the same result in a single pass with a stack and a low-link value \( \mathrm{low}[v] \), the smallest discovery time reachable from \(v\)'s subtree by at most one back or cross edge into a vertex still on the stack. A vertex \(v\) is a component root exactly when \( \mathrm{low}[v] = d[v] \).

Articulation points and bridges

The same low-link machinery finds the vertices and edges whose removal disconnects an undirected graph, which is what network reliability and 2-connectivity questions reduce to. Run DFS and define \( \mathrm{low}[v] = \min \) over \( d[v] \), \( d[w] \) for every back edge \( (v,w) \), and \( \mathrm{low}[c] \) for every child \(c\). Then a non-root vertex \(u\) is an articulation point if and only if it has a child \(v\) with \( \mathrm{low}[v] \ge d[u] \), and the root is one if and only if it has two or more children in the DFS tree. The proof is one observation. Since undirected DFS produces only tree and back edges, the only way \(v\)'s subtree can avoid \(u\) is through a back edge climbing strictly above \(u\), which is precisely \( \mathrm{low}[v] < d[u] \). An edge \( (u,v) \) is a bridge under the strict version, \( \mathrm{low}[v] > d[u] \), because now not even a back edge to \(u\) itself may exist. Both run in \( \Theta(V + E) \), one DFS.

Shortest paths and minimum spanning trees

Dijkstra's algorithm, with the proof and its failure mode

For a directed graph with non-negative edge weights and source \(s\), Dijkstra (1959) maintains a set \(S\) of settled vertices whose distances are final and a tentative distance \( d[v] \) for every vertex, initialized \( d[s] = 0 \), \( d[v] = \infty \). Repeatedly it settles the unsettled vertex \(u\) of minimum \( d[u] \) and relaxes its out-edges, \( d[v] \leftarrow \min(d[v],\, d[u] + w(u,v)) \).

Invariant. When \(u\) is settled, \( d[u] = \delta(s, u) \), the true shortest-path distance. The proof is by induction on settling order. For the base, \(s\) settles with \( d[s] = 0 \), correct since weights are non-negative. For the step, suppose \(u\) is about to settle with tentative \( d[u] \), and suppose toward contradiction some path \(P\) from \(s\) to \(u\) is shorter than \( d[u] \). \(P\) starts inside \(S\) and ends outside, so it has a first edge \( (x, y) \) crossing from \( S \) to \( V \setminus S \). By induction \( d[x] = \delta(s,x) \), and \(y\) was relaxed when \(x\) settled, so

where the middle inequality holds because the prefix of \(P\) through \( (x,y) \) costs at least \( \delta(s,x) + w(x,y) \) and the rest of \(P\), from \(y\) to \(u\), has non-negative weight. But then \( d[y] < d[u] \) with \(y\) unsettled contradicts the choice of \(u\) as the minimum. The italicized clause is the exact point where non-negativity enters, and with one negative edge the argument, and the algorithm, genuinely fail. On edges \( s \to a \) (2), \( s \to b \) (3), \( b \to a \) (\(-2\)), Dijkstra settles \(a\) at distance 2 and never revisits it, while the true distance is \( 3 - 2 = 1 \). Reweighting by potentials repairs negativity when a potential exists. With any \( h : V \to \R \) satisfying \( w(u,v) + h(u) - h(v) \ge 0 \), shortest paths are preserved (the correction telescopes to \( h(s) - h(t) \), the same for every \(s\)-\(t\) path), which is Johnson's algorithm (one Bellman-Ford to find \(h\), then \(n\) Dijkstras) and equally the reason A* with a consistent heuristic is Dijkstra on reweighted edges.

With a binary heap the running time is \( O((V + E) \log V) \), since each vertex is extracted once and each edge causes at most one decrease-key. Fibonacci heaps (Fredman and Tarjan, 1987) make decrease-key \(O(1)\) amortized for \( O(E + V \log V) \), asymptotically better and practically slower except on dense graphs, a constant-factor lesson of the same species as Strassen's.

Bellman-Ford, negative cycles, and the measured gap

Bellman-Ford (Bellman 1958, Ford 1956) handles negative edges by brute democratic relaxation, \( V - 1 \) passes, each relaxing every edge. Claim. After pass \(k\), \( d[v] \) is at most the weight of the shortest path from \(s\) to \(v\) using at most \(k\) edges. The induction on \(k\) is trivially true at \( k = 0 \). A shortest \( \le k \)-edge path to \(v\) is a shortest \( \le k{-}1 \)-edge path to some \(u\) plus edge \( (u,v) \), and pass \(k\) relaxes \( (u,v) \) after \( d[u] \) has reached that prefix value by induction. A shortest path with no negative cycles is simple, hence has at most \( V - 1 \) edges, so \( V - 1 \) passes suffice, for \( O(VE) \) total. A further pass that still relaxes some edge certifies a negative cycle reachable from \(s\), and following parent pointers from the relaxed vertex \(V\) steps back lands inside the cycle. This detection is not a corner case. It is the algorithm's second job, currency arbitrage being the textbook instance (a negative cycle in \( -\log \) exchange rates is a money pump) and distributed distance-vector routing being the historical one.

Measured on random graphs with \( E = 8V \), Dijkstra beat early-exit Bellman-Ford by a stable factor. At \( V = 1600, E = 12{,}800 \), the times were 2.64 ms against 8.80 ms, ratio 3.3, with both algorithms agreeing on every distance on all 100 cross-checked instances (Bellman-Ford was also verified against Floyd-Warshall on 100 instances, zero mismatches). The sharper measurement is about the \( V - 1 \) bound itself. On these random graphs relaxation reached a fixed point after only 6 passes at \( V = 200 \), 8 at 400, and 10 at 800, versus worst-case bounds of 199, 399, 799, so the early-exit variant took 3.31 ms at \( V = 800 \) while the oblivious full \( V-1 \)-pass version took 241.4 ms, a \(73\times\) gap with a fitted exponent of 2.001 for the full version, exactly the \( \Theta(VE) = \Theta(V^2) \) prediction at constant average degree. The worst case is real (a path graph relabeled adversarially needs all passes), but the bound is a guarantee, not a forecast.

Floyd-Warshall: all pairs, by a different subproblem

For all-pairs shortest paths the interesting design choice is what to index the subproblem by. Floyd-Warshall (1962) indexes by the set of permitted intermediate vertices. Let \( d^{(k)}[i][j] \) be the length of the shortest \( i \to j \) path whose interior vertices all lie in \( \{1,\dots,k\} \). Condition on whether vertex \(k\) is used. If it is not, the path is already counted in \( d^{(k-1)}[i][j] \). If it is, then because a shortest path with no negative cycle is simple it uses \(k\) exactly once, and splits into an \( i \to k \) piece and a \( k \to j \) piece, each interior to \( \{1,\dots,k-1\} \) and each optimal by cut-and-paste.

Three lines of loops, \( \Theta(V^3) \) time, \( \Theta(V^2) \) space because the update can be done in place. During round \(k\), row \(k\) and column \(k\) cannot change, since \( d[i][k] \) would have to improve via \(k\) itself, which contributes zero. The \(k\) loop must be outermost, and swapping it inward is the single most common implementation bug in the algorithm. The recurrence, not the code, says why. Negative cycles show up as \( d[i][i] < 0 \), and predecessor matrices reconstruct paths. Against \(V\) runs of Dijkstra at \( O(VE\log V) \), Floyd-Warshall wins on dense graphs and on the transitive-closure variant where min and plus become or and and (Warshall's algorithm), and Johnson's reweighting covers sparse graphs with negative edges. The implementation here agreed with Bellman-Ford on all 100 cross-checked random instances.

Minimum spanning trees: one lemma, three algorithms

Every classical MST algorithm is a corollary of a single exchange lemma. Cut property. For any cut \( (A, V \setminus A) \) and any edge \(e\) that is strictly the lightest crossing the cut, every MST contains \(e\). (Assume distinct weights, with ties handled by any fixed tie-breaking order, which also makes the MST unique.) For the proof, suppose a spanning tree \(T\) omits \(e = (u,v)\). Adding \(e\) to \(T\) closes a cycle, and that cycle crosses the cut an even number of times, in particular at some other edge \( f \ne e \). Then \( T' = T + e - f \) is again a spanning tree (removing an edge of a cycle cannot disconnect), and \( w(T') = w(T) + w(e) - w(f) < w(T) \) since \( w(e) < w(f) \). So \(T\) was not minimum. The companion cycle property (the strictly heaviest edge of any cycle is in no MST) is proved by the same exchange run backward.

Kruskal (1956) sorts edges ascending and adds each edge that joins two different components. The added edge is the lightest edge crossing the cut separating those components from the rest, so the cut property licenses every addition. With union-find under path compression and union by rank the total is \( O(E \log E) \) for the sort plus \( O(E\, \alpha(V)) \) for the merges, the inverse-Ackermann bound from the amortized analysis cited earlier. Prim (1957, also Jarník 1930) grows one component and repeatedly adds the lightest edge leaving it, the cut property applied to the same cut every step. With a binary heap it matches Dijkstra's \( O(E \log V) \) shape. Borůvka (1926), the oldest and the most modern, has every component add its lightest incident edge simultaneously, halving the component count per round for \( O(E \log V) \) total. Its round structure parallelizes naturally, and it is a component of Karger, Klein, and Tarjan's randomized expected-linear-time MST algorithm (1995), while Chazelle's deterministic \( O(E\, \alpha(E,V)) \) algorithm (2000) via soft heaps remains the best known deterministic bound, with the true deterministic complexity of MST still open, one of the cleanest open problems in the field.

Maximum flow and minimum cut

Flows, cuts, and weak duality

A flow network is a directed graph with capacities \( c(u,v) \ge 0 \), a source \(s\), and a sink \(t\). A flow assigns \( f(u,v) \in [0, c(u,v)] \) to each edge with conservation \( \sum_u f(u,v) = \sum_w f(v,w) \) at every internal vertex. Its value \( |f| \) is the net flow out of \(s\). An \(s\)-\(t\) cut is a partition \( (S, T) \) with \( s \in S, t \in T \), of capacity \( c(S,T) = \sum_{u \in S, v \in T} c(u,v) \), counting only forward edges. Weak duality is an accounting identity plus two bounds. Sum conservation over \( v \in S \) to get

dropping the non-negative backward term and capping the forward one. So every flow value is at most every cut capacity, and if some flow's value equals some cut's capacity, both are simultaneously optimal. The content of the max-flow min-cut theorem is that this always happens.

The residual graph and the max-flow min-cut theorem, both directions

The residual graph \( G_f \) encodes every legal local change, edge \( (u,v) \) with residual capacity \( c(u,v) - f(u,v) \) if positive (send more), and reverse edge \( (v,u) \) with residual capacity \( f(u,v) \) if positive (cancel existing flow). An augmenting path is an \( s \to t \) path in \( G_f \). Pushing its bottleneck residual capacity along it yields a valid flow of strictly larger value, cancellation on reverse edges preserving conservation automatically.

Theorem (Ford-Fulkerson 1956). The following are equivalent. (1) \(f\) is a maximum flow. (2) \(G_f\) has no augmenting path. (3) \( |f| = c(S,T) \) for some cut.

(1 \(\Rightarrow\) 2) is the contrapositive of augmentation increasing value. (3 \(\Rightarrow\) 1) is weak duality. The substance is (2 \(\Rightarrow\) 3), where the certificate is constructed. Suppose no augmenting path exists, and let \( S \) be the set of vertices reachable from \(s\) in \( G_f \), so \( t \notin S \) and \( (S, T) \) is a genuine cut. Take any \( u \in S, v \in T \). If \( (u,v) \) is an edge, its residual capacity must be zero, else \(v\) would be reachable, so \( f(u,v) = c(u,v) \), every forward edge saturated. If \( (v,u) \) is an edge, its residual reverse capacity \( f(v,u) \) must be zero, else again \(v\) would be reachable, so every backward edge is empty. Plug both into the accounting identity.

Equality with a cut certifies maximality by weak duality, completing the cycle of implications. The proof is an algorithm. Run any search from \(s\) in \(G_f\) at termination and read off the min cut, which is exactly how the measured verification below extracted cuts. Two corollaries ride along. Integrality. With integer capacities every augmentation increases \( |f| \) by an integer, so a maximum flow exists that is integral, the fact that makes flow encode combinatorics (matchings, disjoint paths) exactly. Termination caveat. With irrational capacities pure Ford-Fulkerson can fail to terminate and even converge to the wrong value. The choice of augmenting path matters, which is the next subsection.

Menger's theorem is the unit-capacity special case worth stating. The maximum number of edge-disjoint \(s\)-\(t\) paths equals the minimum number of edges whose removal disconnects \(t\) from \(s\), immediately from integral max-flow min-cut with all capacities 1.

Edmonds-Karp and the shortest-augmenting-path bound

Edmonds and Karp (1972) fixed the path choice. Always augment along a shortest path in \(G_f\) (BFS). Two lemmas give the bound. Lemma 1. The BFS distance \( d_f(s, v) \) is non-decreasing over augmentations. Augmenting only saturates edges on a shortest path and only creates reverse edges pointing from level \( k+1 \) back to level \( k \). A hypothetical shortening would need a new edge jumping forward more than one level, which augmentation never creates (formally, induction on the first vertex whose distance decreased). Lemma 2. Each edge \( (u,v) \) can be the bottleneck at most \( O(V) \) times. When it bottlenecks, \( d(u) = k \) and \( d(v) = k+1 \), and the edge vanishes from \(G_f\). For it to reappear, flow must be pushed on \( (v,u) \), which requires \( d(v) = k' \) and \( d(u) = k' + 1 \) at that later time, and by Lemma 1 \( k' \ge k + 1 \), so \(u\)'s distance has grown by at least 2 between consecutive bottleneckings. Distances are bounded by \(V\), so each of the \( 2E \) residual edges bottlenecks \( O(V) \) times, each augmentation has at least one bottleneck, and the count of augmentations is \( O(VE) \), each costing one \( O(E) \) BFS.

polynomial with no dependence on capacity values, unlike raw Ford-Fulkerson whose augmentation count can be \( \Theta(|f^*|) \) on the notorious two-triangle example with a capacity-1 cross edge between capacity-\(10^6\) arms, where alternating poorly chosen paths take two million augmentations to move flow two units at a time.

Dinic's algorithm, unit-capacity graphs, and bipartite matching

Dinic's algorithm (1970) batches augmentations. Build the BFS level graph of \(G_f\), find a blocking flow (a flow saturating at least one edge on every \(s\)-\(t\) path of the level graph) with DFS plus retreating pointer bookkeeping, add it, and repeat. Each blocking-flow phase increases \( d_f(s,t) \) by at least 1, so there are at most \(V\) phases, and a phase costs \( O(VE) \) with the pointer trick, giving \( O(V^2 E) \) in general. On matching-type networks, unit capacities with every internal vertex having in-degree or out-degree 1, the analysis sharpens (Even and Tarjan 1975, Karzanov 1973). A phase costs only \( O(E) \) because every saturated edge is dead for the rest of the phase, and after \(k\) phases the residual distance exceeds \(k\), so the remaining flow decomposes into internally vertex-disjoint paths of length \( > k \), of which there can be at most \( V / k \). Each phase adds at least one unit, so at most \( V/k \) phases remain, and choosing \( k = \sqrt{V} \) bounds the total phases by \( 2\sqrt{V} \) for \( O(E \sqrt{V}) \) overall (for general unit edge capacities the same argument with edge-disjoint paths gives \( O(E^{3/2}) \)). Bipartite maximum matching reduces to exactly such a network (source to every left vertex, every right vertex to sink, all capacities 1), where \( O(E \sqrt{V}) \) reproduces Hopcroft-Karp (1973) exactly. Integrality of max flow is what makes the matching integral.

The measurements exercise both the theorem and the algorithm. On 400 random small networks (up to 9 vertices, where exhaustive enumeration of all cuts is feasible), Dinic's flow value equalled the brute-force minimum cut capacity on every instance. On random bipartite graphs, maximum matching via Dinic equalled minimum vertex cover computed independently on all 200 instances, which is Kőnig's theorem (max matching \(=\) min vertex cover in bipartite graphs) observed rather than assumed. Kőnig's theorem itself is max-flow min-cut in a costume, the vertex cover being read off the min cut's saturated structure. On bipartite instances with \( E = 8V_{\text{side}} \), perfect or near-perfect matchings of size 200, 400, 800 (one short, at 799), and 1600 were found in 5.45, 11.61, 34.09, 77.91 ms, a fitted exponent of 1.307 in \(V\) against the \( O(E\sqrt{V}) \sim V^{1.5} \) worst-case shape, the discount coming from random graphs finding augmenting structure faster than adversarial ones.

Push-relabel, and the road to almost-linear time

Goldberg and Tarjan (1988) inverted the augmenting-path viewpoint. Instead of maintaining a feasible flow and pushing toward optimality, push-relabel maintains a preflow (vertices may hold excess) and a height function, pushing excess downhill along admissible edges and relabeling a vertex upward when it is stuck. When no excess remains anywhere, the preflow is a maximum flow. The height function is a certified lower bound on residual distance to \(t\), and the potential-function analysis gives \( O(V^2 E) \) generically and \( O(V^3) \) or \( O(V^2 \sqrt{E}) \) with FIFO or highest-label orderings. With dynamic trees it is \( O(VE \log(V^2/E)) \). For decades this family, plus Goldberg-Rao's \( O(E \min(V^{2/3}, E^{1/2}) \log(V^2/E) \log U) \) bound (1998), was the practical and theoretical frontier, and hand-tuned highest-label push-relabel with gap and global relabeling heuristics remains what serious solvers ship.

The theoretical frontier then moved twice in a decade, driven by a merger of combinatorics with continuous optimization. Interior-point methods recast max flow as a sequence of electrical-flow (Laplacian) solves, each nearly linear after Spielman-Teng. The line through Mądry's work and Kathuria-Liu-Sidford brought unit-capacity max flow to \( E^{4/3+o(1)} \), and then Chen, Kyng, Liu, Peng, Probst Gutenberg, and Sachdeva (2022) achieved maximum flow and minimum-cost flow in \( E^{1+o(1)} \log U \) time, almost linear, via an interior-point method whose updates are served by a dynamic data structure over low-stretch spanning trees. The result, from a team spanning Georgia Tech, ETH Zurich, Stanford, Waterloo, and Toronto, won the FOCS 2022 best paper award, and the follow-up line has pushed toward making the components practical. As of now the almost-linear algorithms are not yet competitive with push-relabel implementations on real instances, a Strassen-shaped gap between the exponent and the constant that anyone who has read this far will find familiar.

Modeling with flow: four reductions

Max flow is valuable less as an algorithm than as a target language. The skill is seeing that a problem with no visible network in it is a flow problem, and the four reductions below are the standard vocabulary. Each has the same shape. Build a network, prove a bijection between the objects being optimized and the cuts or flows of that network, then let a solver do the work.

Bipartite matching and Kőnig's theorem. Orient a bipartite graph \( G = (L \cup R, E) \) as a network, with \( s \to u \) at capacity 1 for each \( u \in L \), \( u \to v \) at capacity \( \infty \) for each edge, and \( v \to t \) at capacity 1 for each \( v \in R \). Integrality makes every max flow a 0/1 assignment, the saturated middle edges form a matching (capacity-1 endpoints force degree at most one on each side), and conversely a matching gives a flow of its size, so max flow \(=\) max matching. Now read the min cut. Let \( (S,T) \) be a minimum cut. Its capacity is finite, so no \( \infty \) edge crosses. Put \( C = (L \cap T) \cup (R \cap S) \), of size exactly \( c(S,T) \) because the crossing edges are precisely \( s \to (L \cap T) \) and \( (R \cap S) \to t \). \(C\) is a vertex cover, because an edge \( (u,v) \) with \( u \in L \cap S \) and \( v \in R \cap T \) would be an \( \infty \) edge crossing the cut. Hence \( \min |{\rm VC}| \le c(S,T) = |f^*| = \) max matching, and the reverse inequality is trivial since a cover needs a distinct vertex per matched edge. Equality is Kőnig's theorem, and the 200-instance check reported above, matching size equal to independently computed minimum vertex cover on every instance, is the theorem observed rather than assumed. The complement gives Kőnig-Egerváry for maximum independent set in bipartite graphs, which is why bipartite instances of an NP-hard problem are easy.

Project selection (maximum-weight closure). Projects have profits and require equipment with costs. A project can only be run if everything it needs is bought, and equipment is shared. Build \( s \to \) project \(i\) with capacity \( p_i \), equipment \( j \to t \) with capacity \( c_j \), and project \( \to \) equipment with capacity \( \infty \). Any finite cut \( (S,T) \) has no \( \infty \) edge crossing, so the selected set \( A = S \cap \text{projects} \) is closed, meaning all equipment it needs is in \(S\). Its capacity is

so minimizing the cut maximizes profit, and the optimal project set is read off the source side of the min cut. For a worked instance, take projects \( P_1 \) worth 100 needing \( M_1 \), \( P_2 \) worth 200 needing \( M_1, M_2 \), and \( P_3 \) worth 150 needing \( M_2, M_3 \), with equipment costs \( M_1 = 100 \), \( M_2 = 100 \), \( M_3 = 200 \), so \( P_{\text{total}} = 450 \). Enumerate the closed sets. \( \{P_1\} \) nets \( 100 - 100 = 0 \), \( \{P_2\} \) nets \( 200 - 200 = 0 \), \( \{P_1, P_2\} \) nets \( 300 - 200 = 100 \) because \( M_1 \) is shared, \( \{P_3\} \) nets \( -150 \), and \( \{P_1,P_2,P_3\} \) nets \( 450 - 400 = 50 \). The optimum is 100. The corresponding cut takes \( S = \{s, P_1, P_2, M_1, M_2\} \) and cuts \( s \to P_3 \) (150), \( M_1 \to t \) (100), and \( M_2 \to t \) (100), for capacity 350, and indeed \( 450 - 350 = 100 \). The identity is exact, not approximate, which is the whole point of a reduction.

Image segmentation. Label each pixel foreground or background with per-pixel likelihoods \( a_i \) (foreground) and \( b_i \) (background) and a separation penalty \( p_{ij} \ge 0 \) for adjacent pixels given different labels. Maximizing \( \sum_{i \in A} a_i + \sum_{j \in B} b_j - \sum_{\text{separated}} p_{ij} \) is the same as minimizing \( \sum_{i \in A} b_i + \sum_{j \in B} a_j + \sum_{\text{separated}} p_{ij} \), since the two objectives sum to the constant \( \sum_i (a_i + b_i) \). That second expression is literally the capacity of the cut in the network with \( s \to i \) of capacity \( a_i \), \( i \to t \) of capacity \( b_i \), and a pair of opposing arcs of capacity \( p_{ij} \) between neighbors. Putting \(i\) on the source side means paying \(b_i\), and separating a neighboring pair means paying \( p_{ij} \) exactly once. This is the Greig-Porteous-Seheult (1989) exact MAP construction for binary Markov random fields, rediscovered for vision as Boykov-Jolly graph cuts, and Kolmogorov and Zabih (2004) characterized which energies are representable this way, exactly the submodular ones, \( E(0,0) + E(1,1) \le E(0,1) + E(1,0) \). Non-submodular terms and more than two labels push the problem back to NP-hard, where alpha-expansion gives a constant-factor guarantee.

Baseball elimination. Team \(z\) is eliminated if it cannot finish first even assuming it wins out and the rest of the season falls its way. Let \( m = w_z + g_z \) be \(z\)'s maximum. Remaining games among the other teams must be distributed, so build a network with a node per remaining pairing \( \{i,j\} \), arc \( s \to \{i,j\} \) of capacity \( g_{ij} \), arcs \( \{i,j\} \to i \) and \( \{i,j\} \to j \) of capacity \( \infty \), and \( i \to t \) of capacity \( m - w_i \), the number of further wins \(i\) can absorb. A feasible schedule keeping everyone at or below \(m\) exists if and only if the max flow saturates every source arc. For a worked instance, take teams \(A\) at 82 wins, \(B\) at 82, \(C\) at 79, and \(D\) at 77 with 6 games left (2 against each of the others), so \( m = 83 \). Remaining among the others are \( g_{AB} = 3 \), \( g_{AC} = 1 \), \( g_{BC} = 1 \), a total of 5 games to place. Sink capacities are \( 83 - 82 = 1 \) for \(A\), 1 for \(B\), \( 83 - 79 = 4 \) for \(C\). The cut \( S = \{s, \{A,B\}, A, B\} \) has capacity \( g_{AC} + g_{BC} + (m - w_A) + (m - w_B) = 1 + 1 + 1 + 1 = 4 < 5 \), so the max flow is at most 4 and \(D\) is eliminated. The min cut is also the human-readable certificate. The subset \( R = \{A, B\} \) has \( (w_A + w_B + g_{AB})/|R| = (82 + 82 + 3)/2 = 83.5 > 83 \), so some team in \(R\) must exceed \(D\)'s ceiling no matter what. That certificate, due to Hoffman and Rivlin, is exactly what the (2 \(\Rightarrow\) 3) direction of the max-flow min-cut proof hands back, and it is the reason to compute a cut rather than only a number. A proof of impossibility is more useful than a failure to find a schedule.

Linear programming, duality, and rounding

Flow is a special case of something larger. A linear program in standard form maximizes \( c^{\mathsf T} x \) subject to \( Ax \le b \), \( x \ge 0 \), and its dual minimizes \( b^{\mathsf T} y \) subject to \( A^{\mathsf T} y \ge c \), \( y \ge 0 \). Weak duality is one line of algebra. For any feasible \(x\) and \(y\),

the first inequality using \( c \le A^{\mathsf T} y \) with \( x \ge 0 \) and the second using \( Ax \le b \) with \( y \ge 0 \). Every dual feasible point is therefore a certificate of optimality-so-far for the primal, which is the same structure as a cut bounding a flow. Strong duality, that the two optima coincide when both are feasible, is von Neumann's and Gale-Kuhn-Tucker's theorem. The proof is a separating-hyperplane argument (Farkas' lemma) and is out of scope here, with a clean treatment in Schrijver's Theory of Linear and Integer Programming. Complementary slackness sharpens it. At optimality \( x_j > 0 \) forces the \(j\)-th dual constraint tight and \( y_i > 0 \) forces the \(i\)-th primal constraint tight, which is how primal-dual approximation algorithms are designed. Max-flow min-cut is exactly LP duality for the flow polytope, and the reason the optimum is integral is that the incidence matrix of a directed graph is totally unimodular, so every basic feasible solution has integer coordinates. Whenever a combinatorial LP has an integral optimum, total unimodularity or a matroid is usually the reason.

When integrality fails, relax and round. Vertex cover as an integer program is \( \min \sum_v x_v \) subject to \( x_u + x_v \ge 1 \) for every edge and \( x_v \in \{0,1\} \). Dropping the integrality constraint to \( x_v \in [0,1] \) gives an LP solvable in polynomial time whose optimum \( \text{OPT}_{\text{LP}} \) is a lower bound on \( \text{OPT} \). Round by setting \( \hat x_v = 1 \) whenever \( x_v \ge 1/2 \). Every edge constraint \( x_u + x_v \ge 1 \) forces at least one endpoint to be at least \(1/2\), so \( \hat x \) is a genuine cover, and \( \sum_v \hat x_v \le \sum_v 2 x_v = 2\,\text{OPT}_{\text{LP}} \le 2\,\text{OPT} \), a second proof of the 2-approximation, this time with a computable lower bound rather than a matching. The gap between the LP and the integer optimum is the integrality gap, and it caps what any LP-rounding argument can achieve. On the complete graph \( K_n \) the all-halves solution has LP value \( n/2 \) while the integer optimum is \( n - 1 \), so the gap approaches 2 and no rounding of this relaxation beats the factor already obtained. Stronger relaxations (semidefinite programs) do better elsewhere, most famously Goemans and Williamson's 0.878-approximation for MAX-CUT, which the Unique Games Conjecture says is optimal.

NP-completeness

The definitions, stated so they can be used

A decision problem is in P if some deterministic algorithm decides it in time polynomial in the input length. It is in NP if yes-instances have certificates verifiable in polynomial time, meaning \( L \in \) NP iff there is a polynomial-time verifier \(V\) and a polynomial \(q\) with \( x \in L \iff \exists w,\, |w| \le q(|x|),\, V(x, w) = 1 \). NP is about checking, not searching. A satisfying assignment, a Hamiltonian cycle, and a factorization are all short and checkable, which is why SAT, Hamiltonicity, and (the decision form of) factoring sit in NP. A polynomial-time reduction \( A \le_p B \) is a polynomial-time map \(f\) with \( x \in A \iff f(x) \in B \). It transfers algorithms backward (\(B\) easy \(\Rightarrow\) \(A\) easy) and hardness forward. \(B\) is NP-hard if every NP problem reduces to it, and NP-complete if additionally \( B \in \) NP. The direction of reduction is the most common working error in the subject. To prove \(B\) hard, reduce a known hard problem to \(B\), never the reverse.

Cook (1971) and, independently, Levin in the Soviet Union proved that SAT is NP-complete. The computation history of any polynomial-time verifier on input \(x\) can be encoded as a polynomial-size Boolean formula (variables for each tape cell, head position, and machine state at each time step, and clauses enforcing legal transitions) that is satisfiable exactly when an accepting certificate exists. The proof is out of scope here beyond that sketch. It lives in Arora and Barak, chapter 2, and in Sipser. Karp (1972) then showed 21 problems complete by reduction chains from SAT, establishing the method that produced the thousands of completeness results catalogued by Garey and Johnson (1979). A reduction in the tradition of that chain is worked in full below.

A complete reduction: 3-SAT to Independent Set

3-SAT asks, given a CNF formula with exactly three literals per clause, whether it is satisfiable. (3-SAT is itself complete, since any clause splits into 3-clauses with fresh variables, \( (a \vee b \vee c \vee d) \equiv (a \vee b \vee z)(\bar z \vee c \vee d) \).) Independent Set asks, given graph \(G\) and integer \(k\), whether \(G\) contains \(k\) pairwise non-adjacent vertices.

Construction. Given \( \varphi \) with clauses \( C_1, \dots, C_m \), build \( G_\varphi \) with one vertex per literal occurrence (\(3m\) vertices), a triangle on the three vertices of each clause, and an edge between every pair of vertices labeled with complementary literals (\(x\) in one clause, \( \bar x \) in another). Set \( k = m \). The construction is computable in time \( O(m^2) \), polynomial.

Satisfiable \( \Rightarrow \) independent set of size \(m\). Fix a satisfying assignment. Each clause has at least one true literal, so pick one such vertex per clause. That is \(m\) vertices. No two share a clause triangle (one per clause), and no two are complementary, because a single assignment cannot make both \(x\) and \( \bar x \) true. So the picked set is independent.

Independent set of size \(m\) \( \Rightarrow \) satisfiable. Let \(I\) be independent with \( |I| = m \). The clause triangles are vertex-disjoint and any independent set contains at most one vertex per triangle, so \(I\) has exactly one vertex in every clause. Define an assignment by setting every literal labeling a vertex of \(I\) to true. This is consistent, because complementary vertices are adjacent and \(I\) is independent, so \(I\) never demands both \(x\) and \( \bar x \). Variables mentioned by no vertex of \(I\) are set arbitrarily. Every clause contains a vertex of \(I\), hence a true literal, so \( \varphi \) is satisfied. Both directions together give \( \varphi \in \text{3-SAT} \iff (G_\varphi, m) \in \text{IS} \), and since IS certificates (the set itself) are checkable in polynomial time, Independent Set is NP-complete.

The chain continues cheaply. For Vertex Cover, \(S\) is a vertex cover iff \( V \setminus S \) is independent (an uncovered edge is exactly an edge inside the complement), so \( (G, k) \in \text{IS} \iff (G, n-k) \in \text{VC} \). For Clique, \(I\) is independent in \(G\) iff \(I\) is a clique in the complement graph \( \bar G \). For Set Cover, given \( (G, k) \) for vertex cover, let the universe be \(E\) and give each vertex \(v\) the set \( S_v = \{ e : v \in e \} \) of edges it touches. A collection of \(k\) sets covering \(E\) is exactly a vertex cover of size \(k\), since covering an edge means picking one of its endpoints. So 3-SAT \( \le_p \) IS \( \le_p \) VC \( \le_p \) Set Cover, a complete chain from a formula to a covering problem, and Set Cover is in NP because a proposed collection is checkable in linear time. Each equivalence is one sentence, which is the point of the method. Hardness spreads through a network of translations, and the practitioner's skill is recognizing that a scheduling conflict graph is asking an independent-set question, or that register allocation is graph coloring (Chaitin's reduction runs both ways), or that a firewall-rule audit is SAT shaped, and then either reaching for a solver or changing the model.

Three standard cautions. NP-hardness is about the worst case of exact solution. It forbids neither good average-case behavior (SAT solvers), nor approximation (next section), nor parameterized tractability (Vertex Cover is solvable in \( O(1.2738^k + kn) \) time, fine for small covers of huge graphs, the fixed-parameter line of Downey and Fellows). Second, NP-completeness is a statement about decision problems under polynomial reductions, not a proof that P \( \ne \) NP. That remains open, and the consensus expectation that they differ rests on decades of failed algorithmic attack plus barrier theorems (relativization, natural proofs, algebrization) explaining why current techniques cannot resolve it. Third, weak versus strong hardness matters operationally. Knapsack's hardness melts under small integer weights, TSP's does not.

Approximation algorithms

Vertex cover: a 2-approximation with a matching lower bound witness

An algorithm is an \( \alpha \)-approximation for a minimization problem if it always outputs a feasible solution of cost at most \( \alpha \cdot \text{OPT} \). The proof obligation is peculiar. OPT is unknown and intractable, so every ratio proof works by exhibiting a computable lower bound on OPT and comparing the algorithm's output to that. For vertex cover, the lower bound is a maximal matching. Greedily pick edges. While any edge has both endpoints uncovered, add the edge to a matching \(M\) and both its endpoints to the cover \(C\). At termination \(C\) is a cover (an uncovered edge would extend \(M\)), and \( |C| = 2|M| \). Any cover must contain at least one endpoint of each edge of \(M\), and edges of \(M\) are disjoint, so \( \text{OPT} \ge |M| \), giving

Measured over 296 random instances (against exact covers from exhaustive search), the mean ratio was 1.6709 and the worst observed was exactly 2.0, attained on instances containing a perfect matching structure. The bound is tight on a single edge already (\(C\) takes 2 vertices, OPT is 1). Whether 2 can be beaten is one of the central questions of hardness of approximation. Dinur and Safra proved \( 1.3606 \)-inapproximability assuming only P \( \ne \) NP, and under Khot's Unique Games Conjecture, Khot and Regev showed \( 2 - \varepsilon \) is hard for every \( \varepsilon \), so the trivial-looking matching algorithm is conjecturally optimal.

Set cover: the greedy \(H_n\) ratio, proved and then measured at its worst

Given a universe \(U\) of \(n\) elements and a family of subsets, cover \(U\) with the fewest sets. Greedy repeatedly takes the set covering the most uncovered elements. Claim. Greedy uses at most \( H_n = 1 + \tfrac12 + \cdots + \tfrac1n \le \ln n + 1 \) times the optimal number \( \text{OPT} \) of sets. The charging argument runs as follows. When greedy picks a set covering \(c\) new elements, charge each of those elements \( 1/c \), so the total charge equals the number of sets greedy uses. Consider the moment an element \(e\) is covered, with \( u \) elements still uncovered just before. The optimal solution covers those \(u\) elements with \( \text{OPT} \) sets, so some set covers at least \( u / \text{OPT} \) of them. Greedy's chosen set covers at least as many, so \(e\)'s charge is at most \( \text{OPT}/u \). Summing over the elements in the order they are covered, the \(j\)-th-from-last element to be covered has \( u \ge j \), so total charge is at most \( \sum_{j=1}^{n} \text{OPT}/j = \text{OPT} \cdot H_n \).

Both halves of this bound were measured. On 271 random instances the mean ratio to the exact optimum was 1.0111 and the worst observed 2.0. On typical instances greedy is nearly optimal, and the logarithm never shows up. It shows up on the adversarial family, which was constructed and run. The universe splits into two halves \(A, B\) (the optimal cover is those 2 sets), while decoy sets \( G_i \) take \( 2^{i-1} \) elements from each half, so greedy, preferring the biggest set at each step (ties broken against it), peels \( G_k, G_{k-1}, \dots \) and uses \( \Theta(\log n) \) sets. In the measurements, universe 16 gives greedy 4 against optimal 2 (ratio 2.0), 64 gives 6 (ratio 3.0), and 512 gives 9 (ratio 4.5), the ratios tracking \( \tfrac{1}{2}\log_2 n = 0.721 \ln n \) against \( \ln n = 6.238 \) at \( n = 512 \). The logarithm is real, and it is also final. Feige (1998) proved that a \( (1 - \varepsilon)\ln n \) approximation for any \( \varepsilon > 0 \) would imply P \(=\) NP-adjacent collapses (quasi-polynomial simulations, with the sharp P \( \ne \) NP form due to Dinur and Steurer, 2014). Greedy set cover is optimal among polynomial algorithms, a complete story of algorithm, analysis, tight instance, and matching hardness.

Metric TSP: tree doubling and Christofides

General TSP admits no constant-factor approximation unless P \(=\) NP (an approximator would detect Hamiltonian cycles by setting non-edges to an enormous weight), so approximation lives in the metric case, where distances satisfy the triangle inequality. Tree doubling computes an MST \(T\), duplicates every edge to make an Eulerian multigraph, takes an Euler tour, and shortcuts repeated vertices. The cost accounting is short. Deleting one edge of the optimal tour leaves a spanning path, hence a spanning tree, so \( w(T) \le \text{OPT} \). The Euler tour costs \( 2 w(T) \), and shortcutting never increases cost, by the triangle inequality applied at each skip. Therefore the tour costs at most \( 2\,\text{OPT} \). Christofides (1976) replaces doubling with a minimum-weight perfect matching on the odd-degree vertices of \(T\) (there are evenly many, by the handshake lemma). The optimal tour restricted to those odd vertices splits into two alternating perfect matchings whose combined cost is at most OPT, so the cheaper one costs at most \( \text{OPT}/2 \), and \( T + M \) is Eulerian with cost at most \( \tfrac32 \text{OPT} \), shortcut as before. Measured over 100 random Euclidean instances against exact Held-Karp optima, tree doubling achieved mean ratio 1.1401 and worst 1.3753, comfortably inside its proved bound of 2, the usual daylight between worst case and typical case. The 3/2 barrier stood for 45 years until Karlin, Klein, and Oveis Gharan (2021, from the University of Washington) beat it by a factor of \( 1.5 - 10^{-36} \) with a randomized algorithm analyzed through the geometry of random spanning trees. The gap to the conjectured 4/3 (the integrality gap of the Held-Karp LP relaxation) remains the field's most famous open ratio.

Parameterized tractability: moving the exponent off \(n\)

A problem is fixed-parameter tractable in a parameter \(k\) if it is solvable in \( f(k) \cdot \text{poly}(n) \) time for some computable \(f\), however horrible. The distinction from \( n^{f(k)} \) is the whole subject. At \( k = 20, n = 10^6 \), \( 2^k n \) is a second, while \( n^{20} \) is never. Vertex cover is the model case, by bounded search tree. Pick any uncovered edge \( (u,v) \). Every cover contains \(u\) or \(v\), so branch on the two choices with \( k \) decremented, giving a binary tree of depth \(k\) and \( O(2^k \cdot E) \) total. Kernelization sharpens it first. Any vertex of degree exceeding \(k\) must be in the cover (otherwise all its neighbors are, exceeding the budget), and after removing those and isolated vertices, a yes-instance has at most \( k^2 \) edges, so the instance is preprocessed down to a size depending only on \(k\). Careful branching rules bring the constant to \( O(1.2738^k + kn) \) (Chen, Kanj, Xia, 2010). Not everything is FPT. Clique and independent set parameterized by solution size are W[1]-hard, the parameterized analogue of NP-hardness in Downey and Fellows' hierarchy, and under the Exponential Time Hypothesis clique genuinely needs \( n^{\Omega(k)} \). Treewidth is the other parameter that pays, since a large family of NP-hard graph problems is linear-time on bounded treewidth by DP over a tree decomposition, which is Courcelle's theorem in its algorithmic form.

Adversary arguments

The decision-tree bound counts leaves. The other way to prove a lower bound is to play the adversary, answering queries so as to keep as many inputs alive as possible and arguing that the algorithm cannot stop until it has asked enough. Finding the maximum of \(n\) elements needs at least \( n - 1 \) comparisons. Call an element a candidate if it has never lost. Every comparison makes at most one element a non-candidate, and the algorithm must reduce \(n\) candidates to one, so it needs \( n - 1 \) comparisons. An adversary who always answers consistently with some surviving assignment forces the count. The refined version is more instructive. Finding both the maximum and the minimum needs \( \lceil 3n/2 \rceil - 2 \) comparisons. The adversary gives each element a state in \( \{ \text{virgin}, \text{won only}, \text{lost only}, \text{both} \} \). The algorithm must end with \( n - 1 \) elements known to have lost at least once and \( n-1 \) known to have won, and only a comparison between two virgins produces two units of that knowledge at once, of which there can be at most \( \lfloor n/2 \rfloor \). Every other comparison yields at most one unit, so the total is at least \( 2(n-1) - \lfloor n/2 \rfloor = \lceil 3n/2 \rceil - 2 \), and the pair-then-compare algorithm attains it. Merging two sorted lists of length \(n\) needs \( 2n - 1 \) comparisons by the same style of argument. On the interleaved input \( a_1 < b_1 < a_2 < b_2 < \cdots \) every adjacent pair must be compared directly, since swapping any uncompared adjacent pair produces a different valid answer the algorithm cannot distinguish. Adversary arguments prove things decision trees cannot, because they exploit the semantics of the queries rather than only their arity.

Fine-grained lower bounds: why quadratic sometimes is the answer

NP-hardness says nothing about problems already in P, where the practical question is whether \( O(n^2) \) can become \( O(n^{1.99}) \). Fine-grained complexity supplies conditional answers by reduction from a small set of conjectures, SETH (the Strong Exponential Time Hypothesis, that CNF-SAT needs \( 2^{(1-o(1))n} \) as clause width grows), the 3SUM conjecture, and the APSP conjecture. Backurs and Indyk (2015, MIT) showed that an \( O(n^{2 - \varepsilon}) \) edit distance algorithm would refute SETH, by encoding a CNF instance's variable-assignment halves into gadget strings whose distance detects a satisfying pair. So the \(O(mn)\) DP above is likely optimal up to subpolynomial factors, and indeed the best known for constant alphabets is \( O(n^2 / \log n) \) (Masek and Paterson, via the Four Russians method). The same program, driven by Vassilevska Williams's group and others, ties longest common subsequence, Fréchet distance, and dynamic-graph problems to the same three pillars. The practical consequence is direction. When a problem is fine-grained-hard, effort goes to approximation (Andoni-Onak-style near-linear approximate edit distance), to special structure, or to parallelism, not to shaving the exponent.

Worked problems

Problems 1 through 5 are distributed through the theory above. These two exercise the reduction and approximation machinery end to end.

Implementation

Every routine below is the code that produced the measurements quoted above, or a direct simplification of it, and each was verified against an independent oracle. Karatsuba and the FFT were checked against CPython's exact big integers, selection against sorting, Dinic against brute-force enumeration of all cuts, and Held-Karp against permutation search. The code printed here was then extracted back out of this page and re-run against those oracles, so the listings and not merely their ancestors are what passed. The tally covers 200 Karatsuba products against CPython's big integers, five FFT products from 1,536 to 24,576 bits, 200 selections against sorting, 400 Dinic flows against exhaustive min-cut enumeration, 200 Held-Karp tours against permutation search, and 200 knapsack, 100 LCS (value and witness), 100 edit distance, and 100 matrix-chain instances against brute force, with zero mismatches anywhere. The C++ listing was compiled and checked on a hand-computed instance. The verification habit is not optional decoration. Divide-and-conquer index arithmetic and residual-graph bookkeeping are exactly the kinds of code where an off-by-one produces plausible wrong answers.

Karatsuba on Python integers, splitting on bit length. The cutoff matters. Below it, the \(\Theta(n^2)\) method's smaller constant wins, as the crossover measurements showed (no speedup until 32 limbs).

def karatsuba(x: int, y: int, cutoff: int = 64) -> int:
    """Multiply non-negative ints; T(n) = 3T(n/2) + O(n) = O(n^1.585)."""
    if x < (1 << cutoff) or y < (1 << cutoff):
        return x * y                     # base case: hardware multiply
    n = max(x.bit_length(), y.bit_length())
    h = n // 2
    mask = (1 << h) - 1
    x1, x0 = x >> h, x & mask            # x = x1 * 2^h + x0
    y1, y0 = y >> h, y & mask
    a = karatsuba(x1, y1, cutoff)        # high * high
    b = karatsuba(x0, y0, cutoff)        # low * low
    m = karatsuba(x1 + x0, y1 + y0, cutoff) - a - b   # the saved product
    return (a << (2 * h)) + (m << h) + b

The iterative radix-2 FFT with bit-reversal permutation, plus polynomial multiplication through the convolution theorem. The recursion of the derivation is unrolled into passes over block sizes 2, 4, 8, ..., which is how every production FFT is written. The recursive version allocates \( \Theta(n \log n) \) temporaries. The rounding step at the end is licensed only by the error analysis. Keep coefficient magnitudes small enough that accumulated error stays far below one half (measured worst case here, \( 4.65 \times 10^{-4} \)).

import cmath

def fft(a, invert=False):
    """In-place iterative Cooley-Tukey; len(a) must be a power of two."""
    n = len(a)
    a = list(a)
    j = 0
    for i in range(1, n):                # bit-reversal permutation
        bit = n >> 1
        while j & bit:
            j ^= bit; bit >>= 1
        j |= bit
        if i < j:
            a[i], a[j] = a[j], a[i]
    length = 2
    while length <= n:                   # butterfly passes
        ang = (2 if not invert else -2) * cmath.pi / length
        wl = cmath.exp(1j * ang)         # primitive length-th root
        for i in range(0, n, length):
            w = 1.0
            for k in range(i, i + length // 2):
                u, v = a[k], a[k + length // 2] * w
                a[k]               = u + v      # e_k + w^k o_k
                a[k + length // 2] = u - v      # e_k - w^k o_k
                w *= wl
        length <<= 1
    if invert:
        a = [x / n for x in a]           # F^{-1} = conj(F)/n
    return a

def poly_mul(p, q):
    """Coefficient product in O(n log n) via the convolution theorem."""
    n = 1
    while n < len(p) + len(q) - 1:
        n <<= 1
    fp = fft(p + [0.0] * (n - len(p)))
    fq = fft(q + [0.0] * (n - len(q)))
    prod = fft([x * y for x, y in zip(fp, fq)], invert=True)
    return [round(c.real) for c in prod[: len(p) + len(q) - 1]]

assert poly_mul([1, 2], [3, 4]) == [3, 10, 8]   # Problem 3, mechanized

Worst-case linear selection. The two recursive calls are visibly on a \( \tfrac15 \)-size and (at most) \( \tfrac{7}{10} \)-size input, the inequality the whole analysis rests on. The base case at 25 elements avoids degenerate group counts.

def select(a, k):
    """k-th smallest (0-indexed), worst-case O(n) by median of medians."""
    if len(a) <= 25:
        return sorted(a)[k]
    groups = [a[i:i + 5] for i in range(0, len(a), 5)]
    medians = [sorted(g)[len(g) // 2] for g in groups]
    pivot = select(medians, len(medians) // 2)     # T(n/5)
    lo = [x for x in a if x < pivot]
    eq = [x for x in a if x == pivot]
    hi = [x for x in a if x > pivot]
    if k < len(lo):                                 # T(7n/10 + 6) at worst
        return select(lo, k)
    if k < len(lo) + len(eq):
        return pivot
    return select(hi, k - len(lo) - len(eq))

Dinic's algorithm. The representation is the load-bearing choice. Arcs are stored in one flat array with arc \( 2e \) and its reverse \( 2e{+}1 \) adjacent, so the residual update after a push is two array writes (cap[id] -= d; cap[id ^ 1] += d) with no hashing and no edge lookup. The it pointers implement the blocking-flow trick. A vertex never rescans arcs it has already exhausted within a phase, which is where the \( O(VE) \) per-phase bound comes from. This is the code that matched brute-force min cuts on all 400 random networks and found the 1600-vertex bipartite matching in 77.91 ms.

class Dinic:
    def __init__(self, n):
        self.n = n
        self.g = [[] for _ in range(n)]   # g[v] = arc ids out of v
        self.to, self.cap = [], []        # arc i and i^1 are paired

    def add_edge(self, u, v, c):
        self.g[u].append(len(self.to)); self.to.append(v); self.cap.append(c)
        self.g[v].append(len(self.to)); self.to.append(u); self.cap.append(0)

    def bfs(self, s, t):                  # build the level graph
        self.level = [-1] * self.n
        self.level[s] = 0
        q = [s]
        for u in q:
            for eid in self.g[u]:
                v = self.to[eid]
                if self.cap[eid] > 0 and self.level[v] < 0:
                    self.level[v] = self.level[u] + 1
                    q.append(v)
        return self.level[t] >= 0

    def dfs(self, u, t, f):               # advance/retreat with pointers
        if u == t:
            return f
        while self.it[u] < len(self.g[u]):
            eid = self.g[u][self.it[u]]
            v = self.to[eid]
            if self.cap[eid] > 0 and self.level[v] == self.level[u] + 1:
                d = self.dfs(v, t, min(f, self.cap[eid]))
                if d > 0:
                    self.cap[eid] -= d
                    self.cap[eid ^ 1] += d     # reverse arc, no lookup
                    return d
            self.it[u] += 1               # arc exhausted for this phase
        return 0

    def max_flow(self, s, t):
        flow = 0
        while self.bfs(s, t):             # at most V phases
            self.it = [0] * self.n
            while True:
                f = self.dfs(s, t, float('inf'))
                if f == 0:
                    break
                flow += f
        return flow

#include <cstdint>
#include <queue>
#include <vector>

// Arc i and i^1 are an edge and its reverse, so a residual update is
// two array writes. int64 capacities: flow values sum, so they overflow
// int32 sooner than people expect.
struct Dinic {
    struct Arc { int to; int64_t cap; };
    std::vector<Arc> arcs;
    std::vector<std::vector<int>> g;     // g[v] = arc ids out of v
    std::vector<int> level, it;
    explicit Dinic(int n) : g(n), level(n), it(n) {}

    void add_edge(int u, int v, int64_t c) {
        g[u].push_back((int)arcs.size()); arcs.push_back({v, c});
        g[v].push_back((int)arcs.size()); arcs.push_back({u, 0});
    }
    bool bfs(int s, int t) {
        std::fill(level.begin(), level.end(), -1);
        std::queue<int> q; q.push(s); level[s] = 0;
        while (!q.empty()) {
            int u = q.front(); q.pop();
            for (int id : g[u])
                if (arcs[id].cap > 0 && level[arcs[id].to] < 0) {
                    level[arcs[id].to] = level[u] + 1;
                    q.push(arcs[id].to);
                }
        }
        return level[t] >= 0;
    }
    int64_t dfs(int u, int t, int64_t f) {
        if (u == t) return f;
        for (int &i = it[u]; i < (int)g[u].size(); ++i) {
            int id = g[u][i], v = arcs[id].to;
            if (arcs[id].cap > 0 && level[v] == level[u] + 1) {
                int64_t d = dfs(v, t, std::min(f, arcs[id].cap));
                if (d > 0) { arcs[id].cap -= d; arcs[id ^ 1].cap += d; return d; }
            }
        }
        return 0;
    }
    int64_t max_flow(int s, int t) {
        int64_t flow = 0, f;
        while (bfs(s, t)) {
            std::fill(it.begin(), it.end(), 0);
            while ((f = dfs(s, t, INT64_MAX)) > 0) flow += f;
        }
        return flow;
    }
};

Held-Karp, the exponential DP from Problem 4 at full scale. The state table is indexed \( [\text{mask}][\text{endpoint}] \). Iterating masks in increasing order is a valid dependency order because a state's predecessors always have smaller masks. The Rust version is the one to reach for at \( n \ge 20 \), where the table has \( 2^{20} \cdot 20 \) entries and Python's interpreter overhead (the measured \( \approx 29 \) ns-scale units per transition) multiplies into minutes.

def held_karp(D):
    """Exact TSP in O(2^n n^2) time, O(2^n n) space."""
    n = len(D)
    FULL = 1 << n
    INF = float('inf')
    g = [[INF] * n for _ in range(FULL)]   # g[S][j]: start 0, visit S, end j
    g[1][0] = 0                            # mask {0}, standing at city 0
    for S in range(1, FULL):
        if not (S & 1):
            continue                       # every tour starts at city 0
        for j in range(n):
            if not (S >> j) & 1 or g[S][j] == INF:
                continue
            for k in range(n):             # extend the path to city k
                if (S >> k) & 1:
                    continue
                nS = S | (1 << k)
                cand = g[S][j] + D[j][k]
                if cand < g[nS][k]:
                    g[nS][k] = cand
    return min(g[FULL - 1][j] + D[j][0] for j in range(1, n))

// Exact TSP in O(2^n n^2) time and O(2^n n) space.
// g[s][j] = cheapest path starting at 0, visiting exactly mask s, ending at j.
fn held_karp(d: &[Vec<u64>]) -> u64 {
    let n = d.len();
    let full = 1usize << n;
    const INF: u64 = u64::MAX / 4;         // headroom so INF + d never wraps
    let mut g = vec![vec![INF; n]; full];
    g[1][0] = 0;                           // mask {0}, standing at city 0
    for s in 1..full {
        if s & 1 == 0 { continue; }        // every tour starts at 0
        for j in 0..n {
            if (s >> j) & 1 == 0 || g[s][j] >= INF { continue; }
            for k in 0..n {
                if (s >> k) & 1 == 1 { continue; }
                let ns = s | (1 << k);
                let cand = g[s][j] + d[j][k];
                if cand < g[ns][k] { g[ns][k] = cand; }
            }
        }
    }
    (1..n).map(|j| g[full - 1][j] + d[j][0]).min().unwrap()
}

Edit distance in two rows of memory, the recurrence from the theory section transcribed. This produced the kitten/sitting table above and matched brute-force edit enumeration on 100 random instances.

def edit_distance(a: str, b: str) -> int:
    """Levenshtein distance in O(mn) time, O(n) space."""
    m, n = len(a), len(b)
    prev = list(range(n + 1))            # D[0][j] = j
    for i in range(1, m + 1):
        cur = [i] + [0] * n              # D[i][0] = i
        for j in range(1, n + 1):
            cost = 0 if a[i - 1] == b[j - 1] else 1
            cur[j] = min(prev[j] + 1,        # delete a[i-1]
                         cur[j - 1] + 1,     # insert b[j-1]
                         prev[j - 1] + cost) # align a[i-1] with b[j-1]
        prev = cur
    return prev[n]

assert edit_distance("kitten", "sitting") == 3
assert edit_distance("intention", "execution") == 5

The rest of the dynamic programs, each a transcription of a recurrence derived above. Three implementation details carry the theory. The knapsack loop runs capacities downward, which is what enforces "each item at most once", since an ascending loop would let \( K[c - w_i] \) already contain item \(i\), silently solving the unbounded variant. The LCS traceback is a second pass over the finished table, not part of the DP, and it terminates because each step strictly decreases \( i + j \). The matrix chain loop iterates over interval lengths rather than endpoints, because \( m[i][j] \) depends only on strictly shorter intervals and no other order makes the dependencies ready.

def knapsack_01(w, v, W):
    """Max value under capacity W. O(nW) time, O(W) space: pseudo-polynomial."""
    dp = [0] * (W + 1)                      # dp[c] after i items = K[i][c]
    for wi, vi in zip(w, v):
        for c in range(W, wi - 1, -1):      # descending: each item used once
            cand = dp[c - wi] + vi
            if cand > dp[c]:
                dp[c] = cand
    return dp[W]

def lcs(a, b):
    """Longest common subsequence and one witness. O(mn) time and space."""
    m, n = len(a), len(b)
    L = [[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]:
                L[i][j] = L[i - 1][j - 1] + 1
            else:
                L[i][j] = max(L[i - 1][j], L[i][j - 1])
    out, i, j = [], m, n                    # traceback: i + j strictly decreases
    while i and j:
        if a[i - 1] == b[j - 1]:
            out.append(a[i - 1]); i -= 1; j -= 1
        elif L[i - 1][j] >= L[i][j - 1]:
            i -= 1
        else:
            j -= 1
    return L[m][n], "".join(reversed(out))

def matrix_chain(p):
    """p holds n+1 dimensions for n matrices; returns min scalar multiplies."""
    n = len(p) - 1
    m = [[0] * n for _ in range(n)]
    for length in range(2, n + 1):          # intervals, shortest first
        for i in range(n - length + 1):
            j = i + length - 1
            m[i][j] = min(m[i][k] + m[k + 1][j] + p[i] * p[k + 1] * p[j + 1]
                          for k in range(i, j))
    return m[0][n - 1]

assert knapsack_01([3, 4, 5], [30, 50, 60], 8) == 90
assert lcs("AGGTAB", "GXTXAYB") == (4, "GTAB")
assert matrix_chain([10, 100, 5, 50]) == 7500      # vs 75000 for the other order

#include <algorithm>
#include <vector>

// All-pairs shortest paths, Theta(V^3) time, Theta(V^2) space.
// The k loop MUST be outermost: after round k, d[i][j] is the best path whose
// interior vertices lie in {0..k}. Swapping the loops is the classic bug.
void floyd_warshall(std::vector<std::vector<long long>> &d) {
    const long long INF = 4e18;
    int n = (int)d.size();
    for (int k = 0; k < n; ++k)
        for (int i = 0; i < n; ++i) {
            if (d[i][k] >= INF) continue;              // no i -> k path yet
            for (int j = 0; j < n; ++j)
                if (d[k][j] < INF && d[i][k] + d[k][j] < d[i][j])
                    d[i][j] = d[i][k] + d[k][j];
        }
    // Afterwards d[i][i] < 0 certifies a negative cycle through vertex i.
}

// 0/1 knapsack, O(nW) time, O(W) space. The descending capacity loop is the
// only thing preventing an item from being taken twice.
long long knapsack(const std::vector<int> &w,
                   const std::vector<long long> &v, int W) {
    std::vector<long long> dp(W + 1, 0);
    for (size_t i = 0; i < w.size(); ++i)
        for (int c = W; c >= w[i]; --c)
            dp[c] = std::max(dp[c], dp[c - w[i]] + v[i]);
    return dp[W];
}

How it is done in practice

What standard libraries actually run

No production system runs textbook mergesort. CPython's list sort is Timsort (Tim Peters, 2002), a mergesort that detects and merges pre-existing runs, adaptive to the partial order real data almost always has. The measured gap on this machine, 27.57 ms against 546.18 ms for pure-Python mergesort at \( n = 160{,}000 \), is mostly C versus interpreter, but Timsort's \( O(n) \) behavior on nearly-sorted input is an asymptotic, not constant, improvement on the workloads that matter. C++ std::sort is introsort (Musser, 1997), quicksort with a depth counter that bails to heapsort at \( 2\log_2 n \) levels, converting quicksort's \(\Theta(n^2)\) tail risk into a worst-case \( O(n \log n) \) guarantee, plus insertion sort below a small threshold, the small-\(n\) constant-factor point made executable. Rust's unstable sort is pdqsort (Orson Peters), which adds pattern detection and deterministic pivot defense. Its stable sort is a Timsort variant. Selection follows the same pattern. std::nth_element is introselect, quickselect with a median-of-medians escape hatch, exactly the hybrid the measured 2× comparison overhead of pure BFPRT predicts.

Big-integer arithmetic is a staircase of the algorithms derived above. CPython switches from schoolbook to Karatsuba around 2,000 bits, which is why its C Karatsuba ran the 24,576-bit product in 2.92 ms. GMP runs the full ladder, schoolbook to Karatsuba to Toom-3 and higher Toom variants to Schönhage-Strassen FFT multiplication, with per-CPU tuned thresholds, and cryptographic libraries stop at the low rungs because RSA-sized operands (2,048 to 4,096 bits) sit exactly in Karatsuba territory, as the crossover table showed. FFT practice belongs to FFTW, whose planner searches over factorization strategies at runtime, and to cuFFT on GPUs. Deep-learning frameworks apply the convolution theorem selectively, using Winograd's minimal-filtering variant (Lavin and Gray, 2016) for small kernels where the FFT's constants lose.

Flow, matching, and routing at scale

Max flow earns its keep in vision and infrastructure. Graph-cut image segmentation solves an \(s\)-\(t\) min cut per image, and the Boykov-Kolmogorov algorithm, an augmenting-path method tuned to grid graphs, beats asymptotically superior push-relabel variants on that workload, another instance of the model-versus-machine theme. Logistics and traffic engineering are min-cost flow. Google's OR-Tools exposes exactly the successive-shortest-path and scaling algorithms the theory prescribes. Road routing does not run Dijkstra raw. Contraction hierarchies (Geisberger et al., 2008) preprocess the graph by contracting vertices in importance order and adding shortcuts, after which continental shortest-path queries take microseconds, a preprocessing/query trade the pure theory only hints at. And the NP-hard modeling layer is a solved engineering problem at surprising scale. Conflict-driven clause learning SAT solvers, and CP-SAT on top of similar machinery, routinely dispatch scheduling instances with millions of clauses, which is why the practical skill hierarchy now starts with recognizing the reduction and writing the model, and only then, if the solver stalls, designing a bespoke algorithm.

The drilling half of this material, the pattern catalog that turns these paradigms into fast interview and contest implementations (two pointers, sliding window, monotonic stacks, union-find, the standard DP templates), lives in the Algorithms section of this site and is not repeated here. The division is deliberate. That section answers "which pattern is this", and this page answers "why is the pattern correct and what does it cost", which are different skills that fail in different ways. An engineer who knows only the patterns writes a correct sliding window and then models a scheduling problem as something NP-hard when a flow formulation existed. An engineer who knows only the theory proves the right bound and then loses an hour to an off-by-one in a residual graph.

The current research frontier

The most consequential shift of the past few years is the merger of combinatorial and continuous methods. The almost-linear max-flow and min-cost-flow algorithm of Chen, Kyng, Liu, Peng, Probst Gutenberg, and Sachdeva (2022) runs an interior-point method whose slowly-changing electrical subproblems are served by dynamic graph data structures. The same toolkit, in work by van den Brand and coauthors, has since been partially derandomized and extended, and the open engineering question is whether any of it can be made competitive with 1988-vintage push-relabel on real instances. The companion surprise came from Bernstein, Nanongkai, and Wulff-Nilsen (2022), negative-weight single-source shortest paths in near-linear time by purely combinatorial means, low-diameter decompositions plus careful use of the Bellman-Ford relaxation this page derived, ending a decades-long gap between Dijkstra's near-linearity and Bellman-Ford's \( O(VE) \).

Matrix multiplication continues its two-track life. The exponent record moved through Alman and Vassilevska Williams (2021) to Alman, Duan, Vassilevska Williams, Xu, Xu, and Zhou (2024), now \( \omega < 2.3716 \), by increasingly refined analysis of the Coppersmith-Winograd tensor, all of it galactic. On the practical track, DeepMind's AlphaTensor (2022) searched tensor decompositions with reinforcement learning and found small-case algorithms including the 47-multiplication \( 4 \times 4 \) result over \( \mathbb{F}_2 \), and AlphaDev (2023) found branch-count improvements to tiny sorting kernels that were merged into LLVM's libc++, machine search operating exactly at the constant-factor layer the asymptotic theory ignores. In approximation, Karlin, Klein, and Oveis Gharan's \( 1.5 - 10^{-36} \) for metric TSP cracked a 45-year barrier and the race toward the conjectured \( 4/3 \) is live. In hardness, the fine-grained program centered on Vassilevska Williams's group at MIT keeps converting folklore beliefs ("edit distance is quadratic") into conditional theorems. A newer thread, algorithms with predictions (Mitzenmacher and Vassilvitskii's survey is the entry point), redesigns classical algorithms to accept possibly-wrong learned advice with provable robustness, learned indexes and predicted-frequency caching being the deployed examples. It is the most direct contact point between this page and modern ML systems.

Open source to read

• networkx/networkx. Readable pure-Python implementations of nearly every algorithm on this page. Start in networkx/algorithms/flow/ and compare the preflow-push and shortest-augmenting-path variants against the derivations here.
• google/or-tools. The industrial modeling layer. ortools/sat/ is CP-SAT, and ortools/graph/ has tight min-cost-flow and assignment implementations.
• scipy/scipy. scipy/sparse/csgraph/ is Dijkstra, Bellman-Ford, and matching at C speed behind a stable API, instructive for the layer between textbook and production.
• kth-competitive-programming/kactl. The KTH team reference, with tersely correct implementations, each carrying a stated complexity and test. content/graph/Dinic.h is 40 lines.
• cp-algorithms/cp-algorithms. Derivation-first articles with code for FFT, flows, and DP optimizations (divide-and-conquer DP, Knuth's optimization) that go beyond this page.
• atcoder/ac-library. Battle-tested C++ for max flow, min-cost flow, and union-find. atcoder/maxflow.hpp is the paired-arc representation used above.
• boostorg/graph. The generic-programming take. push_relabel_max_flow.hpp shows what the Goldberg-Tarjan heuristics (gap, global relabeling) look like in earnest.
• Z3Prover/z3. The SMT solver that makes "reduce it to SAT and stop writing algorithms" a real strategy. Open src/sat/sat_solver.cpp to see conflict-driven clause learning, the machinery behind the practical collapse of NP-hardness.
• TheAlgorithms/C-Plus-Plus. A broad, uneven, well-commented catalog. graph/dinic_maxflow.cpp is a useful second opinion to diff against the version above when debugging arc bookkeeping.
• orlp/pdqsort. Pattern-defeating quicksort, the basis of Rust's unstable sort. The README is a thorough study in adversarial input analysis at the constant-factor layer.

Common misconceptions

"Big-O tells you which algorithm is faster." It bounds growth rates, one-sidedly, in a stated machine model, as \( n \to \infty \). It is silent at any particular \(n\), silent about constants, and silent about memory hierarchies. The measurements above show Karatsuba losing to schoolbook until 32 limbs, Strassen's wall-clock advantage evaporating while its multiplication count tracks \( n^{2.807} \) perfectly, and Timsort beating an identically-asymptotic mergesort by \(20\times\). The correct reading of \( f = O(g) \) is set membership under quantifiers, nothing more.

"Randomized algorithms work because inputs are usually random." Backwards. Average-case analysis assumes random inputs. Randomized algorithms make no input assumption at all and instead randomize their own decisions so that no fixed input is bad. The adversary, including the denial-of-service attacker feeding sorted arrays to quicksort or colliding keys to a hash table, cannot see the coins. That distinction is why the quicksort expectation \( 2(n+1)H_n - 4n \) holds for every input, and why it matched measurement within 2%.

"Greedy is a heuristic, and if the local choice looks sensible, it is probably near-optimal." Greedy algorithms are either provably exact (interval scheduling, Huffman, MST, anything matroidal), provably approximate with a known ratio (set cover's \( H_n \), which is real and was measured growing), or wrong by an unbounded factor. Three of the four natural greedy rules for interval scheduling fail on five-interval counterexamples. The proof, stays-ahead or exchange, is the deliverable, not the rule.

"Dynamic programming is just caching." Memoization is the implementation. The algorithm is the theorem that the problem decomposes over a polynomial family of subproblems, proved by cut-and-paste. Where that proof fails, as for longest simple path, caching computes confidently wrong answers. The design work is choosing state so that the past interacts with the future only through it, which is why Held-Karp's state is a set and an endpoint, not an ordering.

"NP-complete means give up." NP-hardness constrains worst-case exact polynomial algorithms and nothing else. Vertex cover is FPT in the cover size, knapsack has an FPTAS, set cover greedy is \( \ln n \)-optimal, and CP-SAT closes million-variable scheduling instances daily. The practical meaning of a completeness proof is a licensing decision, to stop searching for an exact polynomial algorithm and pick among solvers, parameters, and ratios.

"Bellman-Ford takes \(V - 1\) passes." That is the guarantee, not the behavior. On the measured random graphs relaxation converged in 6 to 10 passes where the bound said 199 to 799, and the early-exit check turned a 241 ms run into 3.3 ms. Worst-case bounds are ceilings that adversarial instances attain (a badly-ordered path graph does need all passes), not forecasts.

"Max-flow is a networking tool." Max-flow min-cut is a combinatorial duality theorem, the integral case of LP duality, and its applications are mostly not about pipes. Bipartite matching, Kőnig's theorem, Menger's disjoint paths, image segmentation, project selection with prerequisites, and baseball elimination are all min-cut statements. The measured Kőnig verification (200 graphs, zero mismatches) is duality observed empirically.

"Asymptotically optimal algorithms are what production systems run." Fibonacci heaps, BFPRT with groups of five, Strassen at practical sizes, the \( \omega < 2.3716 \) matrix multiplication line, and the almost-linear max-flow algorithms are all, today, outrun by asymptotically worse competitors with better constants and cache behavior. The mature position holds both truths. Exponents decide eventually, constants decide today, and the crossover point is an empirical question the tables above answer for six different algorithms.

Self-check

References

1. Cormen, Leiserson, Rivest, Stein. Introduction to Algorithms, 4th ed. MIT Press, 2022.
2. Kleinberg, Tardos. Algorithm Design. Addison-Wesley, 2006.
3. Dasgupta, Papadimitriou, Vazirani. Algorithms. McGraw-Hill, 2008.
4. Williamson, Shmoys. The Design of Approximation Algorithms. Cambridge University Press, 2011. designofapproxalgs.com
5. Erickson, J. Algorithms. Self-published, 2019. jeffe.cs.illinois.edu/algorithms
6. Vazirani, V. V. Approximation Algorithms. Springer, 2001.
7. Motwani, Raghavan. Randomized Algorithms. Cambridge University Press, 1995.
8. Arora, Barak. Computational Complexity: A Modern Approach. Cambridge University Press, 2009.
9. Garey, Johnson. Computers and Intractability: A Guide to the Theory of NP-Completeness. W. H. Freeman, 1979.
10. Karatsuba, Ofman. Multiplication of multidigit numbers on automata. Doklady Akademii Nauk SSSR 145, 1962 (English translation in Soviet Physics Doklady 7, 1963).
11. Cooley, Tukey. An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation 19, 1965.
12. Strassen. Gaussian elimination is not optimal. Numerische Mathematik 13, 1969.
13. Blum, Floyd, Pratt, Rivest, Tarjan. Time bounds for selection. Journal of Computer and System Sciences 7(4), 1973.
14. Huffman. A method for the construction of minimum-redundancy codes. Proceedings of the IRE 40(9), 1952. doi:10.1109/JRPROC.1952.273898
15. Dijkstra. A note on two problems in connexion with graphs. Numerische Mathematik 1, 1959. doi:10.1007/BF01386390
16. Christofides. Worst-case analysis of a new heuristic for the travelling salesman problem. Report 388, Carnegie Mellon University, 1976, reprinted in Operations Research Forum 3, 2022.
17. Ford, Fulkerson. Maximal flow through a network. Canadian Journal of Mathematics 8, 1956.
18. Edmonds, Karp. Theoretical improvements in algorithmic efficiency for network flow problems. Journal of the ACM 19(2), 1972.
19. Dinic. Algorithm for solution of a problem of maximum flow in a network with power estimation. Soviet Mathematics Doklady 11, 1970.
20. Goldberg, Tarjan. A new approach to the maximum-flow problem. Journal of the ACM 35(4), 1988.
21. Chen, Kyng, Liu, Peng, Probst Gutenberg, Sachdeva. Maximum flow and minimum-cost flow in almost-linear time. FOCS 2022. arXiv:2203.00671
22. Cook. The complexity of theorem-proving procedures. STOC 1971.
23. Karp. Reducibility among combinatorial problems. In Complexity of Computer Computations, Plenum Press, 1972.
24. Karger. Global min-cuts in RNC, and other ramifications of a simple min-cut algorithm. SODA 1993.
25. Tarjan. Efficiency of a good but not linear set union algorithm. Journal of the ACM 22(2), 1975.
26. Feige. A threshold of ln n for approximating set cover. Journal of the ACM 45(4), 1998.
27. Backurs, Indyk. Edit distance cannot be computed in strongly subquadratic time (unless SETH is false). STOC 2015. arXiv:1412.0348
28. Karlin, Klein, Oveis Gharan. A (slightly) improved approximation algorithm for metric TSP. STOC 2021. arXiv:2007.01409
29. Harvey, van der Hoeven. Integer multiplication in time O(n log n). Annals of Mathematics 193(2), 2021.
30. Fawzi, Balog, Huang, Hubert, Romera-Paredes, et al. Discovering faster matrix multiplication algorithms with reinforcement learning. Nature 610, 2022. nature.com/articles/s41586-022-05172-4
31. Bernstein, Nanongkai, Wulff-Nilsen. Negative-weight single-source shortest paths in near-linear time. FOCS 2022. arXiv:2203.03456

Key takeaway