The Engineer's Compass

Propositions and connectives

A proposition is a statement that is either true or false: "7 is prime", "the build passed", "x > 3" once x has a value. Propositions combine with a few connectives, each defined by a truth table: not ($\neg p$), and ($p\wedge q$), or ($p\vee q$), exclusive or ($p\oplus q$), implies ($p\to q$) and if and only if ($p\leftrightarrow q$). Most map straight onto code: not, and, or, and != on booleans for exclusive or. Mathematical "or" is always inclusive. "p or q" is true when both are true, unlike the "or" of a restaurant menu.

Two formulas are equivalent if they have the same truth value in every row of the truth table, that is, for every assignment of true and false to their variables. A formula that is true in every row is valid, or a tautology. One that is true in at least one row is satisfiable. With n variables there are 2ⁿ rows, so a truth table settles any question about a few variables and becomes useless at a few dozen. Deciding satisfiability in general is the SAT problem, the first problem proved NP-complete, which The Engineer's Workbench covers in its last chapter.

Implication, the connective everyone gets wrong

$p\to q$, read "if p then q" or "p implies q", is false in exactly one case: p is true and q is false. In every other row it's true, including both rows where p is false. That surprises people, but it's the only definition that works. "If a number is divisible by 4, it's even" should be true as a statement about all numbers, so it has to be true for 3 (false premise, false conclusion) and for 6 (false premise, true conclusion) as well as for 8. An implication with a false premise is vacuously true. Python follows the same rule: all([]) is True, because a claim about every item of an empty list has no item that could break it.

English hides implications in several phrasings. "p only if q" means $p\to q$: q is necessary for p. "p if q" means $q\to p$: q is sufficient for p. "p unless q" means $\neg q\to p$. Every implication has three relatives, made by swapping its sides, negating both, or both at once.

The contrapositive is the most useful rewriting in mathematics, because proving $\neg q\to \neg p$ is often easier than proving $p\to q$ directly (chapter 2). The converse is the most common error in reasoning: from "if the server is overloaded, requests time out" you can't conclude "requests are timing out, so the server is overloaded".

Checking equivalences by machine

With a few variables you never need to argue about an equivalence. Enumerate every row and compare. The example does that for the four conditionals, De Morgan's laws, and two rewritings of nested implications, then evaluates some quantified statements over a small domain, where "for all" is a loop with all and "there exists" is a loop with any.

from itertools import product

def implies(p, q): return (not p) or q

def check(name, f, g, n):
    """Compare two formulas on every row of the truth table; report the first row where they differ."""
    for row in product((False, True), repeat=n):
        if f(*row) != g(*row):
            vals = ", ".join(f"{v}={int(b)}" for v, b in zip("pqr", row))
            return print(f"   NOT equivalent  {name}   (differ at {vals})")
    print(f"   equivalent      {name}   ({2 ** n} rows checked)")

print("p q | p -> q")
for p, q in product((False, True), repeat=2):
    print(f"{int(p)} {int(q)} |   {int(implies(p, q))}")

print("\nclaimed equivalences:")
check("p -> q  vs  not q -> not p   (contrapositive)", implies, lambda p, q: implies(not q, not p), 2)
check("p -> q  vs  q -> p           (converse)", implies, lambda p, q: implies(q, p), 2)
check("p -> q  vs  not p -> not q   (inverse)", implies, lambda p, q: implies(not p, not q), 2)
check("q -> p  vs  not p -> not q   (converse vs inverse)", lambda p, q: implies(q, p), lambda p, q: implies(not p, not q), 2)
check("p -> q  vs  not p or q", implies, lambda p, q: (not p) or q, 2)
check("not (p and q)  vs  not p or not q   (De Morgan)", lambda p, q: not (p and q), lambda p, q: (not p) or (not q), 2)
check("not (p -> q)  vs  p and not q", lambda p, q: not implies(p, q), lambda p, q: p and not q, 2)
check("p -> (q -> r)  vs  (p and q) -> r", lambda p, q, r: implies(p, implies(q, r)), lambda p, q, r: implies(p and q, r), 3)
check("(p -> q) -> r  vs  p -> (q -> r)", lambda p, q, r: implies(implies(p, q), r), lambda p, q, r: implies(p, implies(q, r)), 3)

D = range(10)                                   # quantifiers over a finite domain are loops
print("\nover the domain 0..9:")
print("   for all x, exists y: x + y = 9   ", all(any(x + y == 9 for y in D) for x in D))
print("   exists y, for all x: x + y = 9   ", any(all(x + y == 9 for x in D) for y in D))
print("   for all x, exists y: x < y       ", all(any(x < y for y in D) for x in D), "(true over all integers)")
P = lambda x, y: x * y % 7 == 1
lhs = not all(any(P(x, y) for y in D) for x in D)
rhs = any(all(not P(x, y) for y in D) for x in D)
print(f"   not (for all x, exists y: xy mod 7 = 1) is {lhs}; exists x, for all y: xy mod 7 != 1 is {rhs}")
print("   all([]) =", all([]), "  any([]) =", any([]), "  (a claim about no items is vacuously true)")

p q | p -> q
0 0 |   1
0 1 |   1
1 0 |   0
1 1 |   1

claimed equivalences:
   equivalent      p -> q  vs  not q -> not p   (contrapositive)   (4 rows checked)
   NOT equivalent  p -> q  vs  q -> p           (converse)   (differ at p=0, q=1)
   NOT equivalent  p -> q  vs  not p -> not q   (inverse)   (differ at p=0, q=1)
   equivalent      q -> p  vs  not p -> not q   (converse vs inverse)   (4 rows checked)
   equivalent      p -> q  vs  not p or q   (4 rows checked)
   equivalent      not (p and q)  vs  not p or not q   (De Morgan)   (4 rows checked)
   equivalent      not (p -> q)  vs  p and not q   (4 rows checked)
   equivalent      p -> (q -> r)  vs  (p and q) -> r   (8 rows checked)
   NOT equivalent  (p -> q) -> r  vs  p -> (q -> r)   (differ at p=0, q=0, r=0)

over the domain 0..9:
   for all x, exists y: x + y = 9    True
   exists y, for all x: x + y = 9    False
   for all x, exists y: x < y        False (true over all integers)
   not (for all x, exists y: xy mod 7 = 1) is True; exists x, for all y: xy mod 7 != 1 is True
   all([]) = True   any([]) = False   (a claim about no items is vacuously true)

The contrapositive matches in all four rows. The converse and the inverse each fail at p = 0, q = 1, where the implication is vacuously true and they're false, but they match each other, as the figure says. The negation of an implication, $\neg (p\to q)\equiv p\wedge \neg q$, is worth memorising: an implication fails only when its premise holds and its conclusion doesn't, which is exactly what a counterexample is. Implication isn't associative. $(p\to q)\to r$ and $p\to (q\to r)$ differ when everything is false. The second form equals $(p\wedge q)\to r$, which is why a chain of conditions in code can be written either as nested ifs or as one if with and.

Predicates and quantifiers

A predicate is a statement with variables, such as P(x, y): "x + y = 9". It becomes a proposition when you fix the variables or quantify them: $\forall x$, "for all x", or $\exists y$, "there exists a y". Quantifiers always range over a domain, and the domain changes the truth value. "For every x there's a larger y" is true over the integers and false over 0..9, where nothing is larger than 9, as the example's third line shows.

The order of mixed quantifiers matters. $\forall x\exists y:x+y=9$ says every x has its own partner, which is true. $\exists y\forall x:x+y=9$ says one y works for every x at once, which is false. Swapping $\forall$ and $\exists$ changes the claim. Specifications of systems are full of this distinction: "every request is served by some replica" is not "some replica serves every request".

To negate a quantified statement, push the negation inward and flip each quantifier as it passes:

This is De Morgan's law for quantifiers, since "for all" is a big "and" and "there exists" is a big "or". The example checks one instance: it's false that every x in 0..9 has a y with xy mod 7 = 1, and correspondingly there's an x, namely 0 (or 7), for which every y fails. Negation is how you find out what would disprove a claim. The negation of the service-level objective "every request gets a response within one second" is "some request doesn't get a response within one second", and that single request is what your monitoring has to catch.

What a proof is

A proof starts from definitions and facts already established and reaches a conclusion by steps that each follow by logic alone. Definitions do most of the work, so they have to be exact. An integer n is even if n = 2k for some integer k, and odd if n = 2k + 1 for some integer k. A real number is rational if it equals p/q for integers p and q with q ≠ 0. Almost every first proof consists of unpacking a definition, doing algebra, and packing the result back into a definition.

Direct proof

To prove $P\to Q$ directly, assume P and derive Q. Claim: if n is odd, then n² is odd. Proof: since n is odd, n = 2k + 1 for some integer k. Then

which is twice an integer plus one, so n² is odd. ∎ That's the whole pattern: definition in, algebra, definition out. The square at the end marks where the proof stops.

Contrapositive and contradiction

Claim: if n² is even, then n is even. A direct proof stalls: knowing n² = 2k says nothing obvious about n. The contrapositive, "if n is odd, then n² is odd", is exactly what we just proved, and chapter 1 showed that an implication and its contrapositive are equivalent. Done.

Proof by contradiction assumes the statement is false and derives something impossible. The classic example is that $\sqrt{2}$ is irrational. Suppose it were rational, $\sqrt{2}=p/q$ with the fraction in lowest terms. Squaring gives $p^2=2q^2$, so p² is even, so p is even by the claim above: p = 2m. Then $4m^2=2q^2$, so $q^2=2m^2$ and q is even too. Both even contradicts lowest terms, so no such fraction exists. ∎

Euclid's proof that there are infinitely many primes has the same shape. Suppose $p_1,p_2,\dots ,p_k$ were all the primes, and let $N=p_1{p}_2\cdots p_k+1$. Dividing N by any $p_i$ leaves remainder 1, so no $p_i$ divides N. But N > 1 has some prime factor (chapter 3 proves that every integer above 1 does), and it isn't on the list. So the list wasn't complete. ∎ Note that N itself needn't be prime. 2 × 3 × 5 × 7 × 11 × 13 + 1 = 30,031 = 59 × 509. The proof only needs N to have a prime factor that's missing from the list.

Cases and existence

Proof by cases splits the situation into possibilities that together cover everything, and proves the claim in each. To show that n² + n is always even: if n is even, both terms are even; if n is odd, both terms are odd and their sum is even. Every integer is even or odd, so the claim holds. The danger in a case split is a missing case, which is why "without loss of generality" deserves a second look whenever you see it.

To prove that something exists, the most convincing way is to construct it. A non-constructive proof shows existence without producing the object. Here's a famous one: there are irrational numbers a and b with a^b rational. Consider ${\sqrt{2}}^{\sqrt{2}}$. If it's rational, take $a=b=\sqrt{2}$. If it's irrational, take $a={\sqrt{2}}^{\sqrt{2}}$ and $b=\sqrt{2}$, since $a^b={\sqrt{2}}^2=2$. Either way such a pair exists, and the proof never tells you which case holds. To disprove a universal claim, one counterexample is enough, and finding counterexamples is where a computer earns its keep.

Examples aren't proofs

The example tests five claims that each held for every case their first observers checked. It uses nothing more sophisticated than trial division.

from math import isqrt

def smallest_factor(n):
    if n % 2 == 0: return 2
    for d in range(3, isqrt(n) + 1, 2):
        if n % d == 0: return d
    return n
def is_prime(n): return n >= 2 and smallest_factor(n) == n

# 1. Euler's prime-generating polynomial
n = next(n for n in range(10**6) if not is_prime(n * n + n + 41))
print(f"n^2 + n + 41 is prime for n = 0 to {n - 1}, then n = {n} gives {n*n + n + 41:,} = 41 x 41")

# 2. Fermat's conjecture that every 2^(2^k) + 1 is prime
for k in range(6):
    F = 2 ** 2 ** k + 1
    d = smallest_factor(F)
    print(f"F{k} = {F:>13,}  " + ("prime" if d == F else f"= {d} x {F // d:,}"))

# 3. "Goldbach's other conjecture": every odd composite is a prime plus twice a square
fails = [n for n in range(9, 100_000, 2) if not is_prime(n)
         and not any(is_prime(n - 2 * k * k) for k in range(1, isqrt(n // 2) + 1))]
print(f"odd composites below 100,000 that are not a prime plus twice a square: {fails}")

# 4. 2^p - 1 for prime p
p = next(p for p in range(2, 100) if is_prime(p) and not is_prime(2 ** p - 1))
print(f"2^p - 1 is prime for p = 2, 3, 5, 7, but 2^{p} - 1 = {2**p - 1:,} = {smallest_factor(2**p - 1)} x {(2**p - 1) // smallest_factor(2**p - 1)}")

# 5. Evidence is not proof, but it can agree with one: no fraction p/q squares to exactly 2
best = []
p, q = 1, 1
while q < 10**6:
    best.append((p, q, p * p - 2 * q * q))
    p, q = p + 2 * q, p + q                     # each step is the next best approximation to sqrt(2)
print("best fractions for sqrt(2), with p^2 - 2q^2:", ", ".join(f"{p}/{q} ({d:+d})" for p, q, d in best[-4:]))

n^2 + n + 41 is prime for n = 0 to 39, then n = 40 gives 1,681 = 41 x 41
F0 =             3  prime
F1 =             5  prime
F2 =            17  prime
F3 =           257  prime
F4 =        65,537  prime
F5 = 4,294,967,297  = 641 x 6,700,417
odd composites below 100,000 that are not a prime plus twice a square: [5777, 5993]
2^p - 1 is prime for p = 2, 3, 5, 7, but 2^11 - 1 = 2,047 = 23 x 89
best fractions for sqrt(2), with p^2 - 2q^2: 47321/33461 (-1), 114243/80782 (+1), 275807/195025 (-1), 665857/470832 (+1)

Euler's polynomial n² + n + 41 produces primes for the forty values 0 to 39, then fails at n = 40, where it gives 41². Fermat believed every number of the form 2^2k + 1 was prime. The first five are, and in 1732 Euler found that the sixth, 4,294,967,297, is 641 × 6,700,417. No Fermat number beyond the fifth has ever been found to be prime. Goldbach suggested to Euler that every odd composite number is a prime plus twice a square, which holds for every odd composite up to 5,775 and then fails at 5,777 and 5,993. Below 100,000 the example finds no other failures, which is also not a proof that there are none. And 2^p − 1 for prime p, the form of every Mersenne prime, is prime for the first four primes and fails at p = 11.

The last line shows the opposite situation, a search that agrees with a proof. The best fractions for $\sqrt{2}$ come ever closer, but $p^2-2q^2$ alternates between $-1$ and $+1$ and never reaches 0, as it would have to for $p/q=\sqrt{2}$. The proof above explains why it never will. The search alone, however long it ran, could never say so.

Programs are the same. Testing a function on a thousand inputs shows that it works on those thousand; a proof of a loop invariant (chapter 3) shows that it works on all of them. Tests are how you find the counterexample when the "proof" in your head was wrong, so you need both.

The principle

Let P(n) be a statement about a natural number n. If P(0) is true (the base case), and P(n) implies P(n + 1) for every n (the inductive step), then P(n) is true for every n. Each case knocks over the next, like dominoes. The assumption P(n) inside the step is the inductive hypothesis, and the whole skill lies in using it.

Claim: for every n ≥ 1,

Proof. Base case: for n = 1 both sides equal 1. Inductive step: assume the formula holds for n. Then

which is the formula for n + 1. ∎ Notice what induction does and doesn't give you. It verifies a formula you already have, but it doesn't find one. Guessing comes first, often from computing small cases, and chapter 11 has methods for finding closed forms directly.

Induction is recursion

A recursive function and an induction proof have the same skeleton. The base case of the recursion is the base case of the proof. The recursive call on a smaller input is the inductive hypothesis, and the work done after the call returns is the inductive step. To prove a recursive function correct, you prove "it returns the right answer for every input of size n" by induction on n, assuming the recursive calls return the right answers because they're on smaller inputs. If you can't write that proof, the function usually has a bug, most often a missing or wrong base case.

Strong induction and well-ordering

Sometimes P(n + 1) depends on cases further back than P(n). Strong induction lets the step assume P(k) for every k up to n. It proves no more than ordinary induction, but it's often far more convenient. Claim: every integer n ≥ 2 is a product of primes. Proof. If n is prime, it's a product of one prime. Otherwise n = ab with 2 ≤ a, b < n, and by the strong hypothesis both a and b are products of primes, so n is too. ∎ This is the step chapter 2's proof that there are infinitely many primes relied on, and the example runs it as a recursive function.

A variant: every amount of postage from 12 cents up can be made from 4-cent and 5-cent stamps. Check the four base cases 12, 13, 14 and 15 directly. For any larger amount n, the amount n − 4 is at least 12, so by the strong hypothesis it can be made, and one more 4-cent stamp makes n. The step reaches back four cases, which is why it needs four base cases.

An equivalent principle is well-ordering: every nonempty set of natural numbers has a smallest element. It gives proofs by smallest counterexample. If the claim failed somewhere, there would be a smallest n where it fails; you then show that n's failure forces a failure at some smaller number, which is impossible.

Induction, run

The example checks two sum formulas on a thousand cases, runs the factorisation and postage proofs as recursive functions, asserts a loop invariant on every iteration of fast exponentiation, and then explores every reachable position of the 8-puzzle to test an invariant of a state machine.

from collections import deque
from math import factorial

# 1. Checking a formula on many cases is evidence; the proof by induction covers the rest.
ok = all(sum(range(1, n + 1)) == n * (n + 1) // 2 for n in range(1, 1001))
ok2 = all(sum(i * i for i in range(1, n + 1)) == n * (n + 1) * (2 * n + 1) // 6 for n in range(1, 1001))
print(f"1 + ... + n = n(n+1)/2 for n = 1..1000: {ok};  sum of squares = n(n+1)(2n+1)/6: {ok2}")

# 2. Strong induction, run as recursion: every n >= 2 is a product of primes.
calls = 0
def factor(n):
    global calls; calls += 1
    for d in range(2, int(n ** 0.5) + 1):
        if n % d == 0:
            return factor(d) + factor(n // d)    # both parts are smaller: the strong hypothesis
    return [n]                                 # no divisor: n itself is prime
print(f"factor(360) = {factor(360)} in {calls} calls; factor(9973) = {factor(9973)}")

def stamps(n):                                 # every amount from 12 cents up, with 4c and 5c stamps
    base = {12: (3, 0), 13: (2, 1), 14: (1, 2), 15: (0, 3)}
    if n in base: return base[n]
    fours, fives = stamps(n - 4)               # strong induction: the case 4 below, plus one stamp
    return fours + 1, fives
print("postage:", ", ".join(f"{n}c = {a}x4 + {b}x5" for n in (16, 23, 31, 100) for a, b in [stamps(n)]))

# 3. A loop invariant: result * base**exp == x**n holds before and after every iteration.
def power(x, n):
    result, base, exp, mults = 1, x, n, 0
    while exp > 0:
        assert result * base ** exp == x ** n
        if exp & 1: result *= base; mults += 1
        base *= base; exp >>= 1; mults += 1
    assert result * base ** exp == x ** n and exp == 0   # invariant + exit condition = correctness
    return result, mults
r, m = power(3, 1000)
print(f"3^1000 by repeated squaring: {m} multiplications instead of 999, invariant held at every step, "
      f"correct: {r == 3 ** 1000}")

# 4. A state-machine invariant: the 8-puzzle. Every slide preserves the parity of inversions.
def inversions(state):
    tiles = [t for t in state if t != 0]
    return sum(tiles[i] > tiles[j] for i in range(len(tiles)) for j in range(i + 1, len(tiles)))
goal = (1, 2, 3, 4, 5, 6, 7, 8, 0)
seen, queue = {goal}, deque([goal])
while queue:
    s = queue.popleft(); z = s.index(0); r, c = divmod(z, 3)
    for dr, dc in ((1, 0), (-1, 0), (0, 1), (0, -1)):
        if 0 <= r + dr < 3 and 0 <= c + dc < 3:
            t = list(s); w = (r + dr) * 3 + c + dc; t[z], t[w] = t[w], t[z]; t = tuple(t)
            if t not in seen: seen.add(t); queue.append(t)
print(f"8-puzzle: {len(seen):,} of {factorial(9):,} arrangements reachable from the goal; "
      f"all with an even number of inversions: {all(inversions(s) % 2 == 0 for s in seen)}")
print(f"swapping tiles 7 and 8 gives {inversions((1, 2, 3, 4, 5, 6, 8, 7, 0))} inversion: reachable: "
      f"{(1, 2, 3, 4, 5, 6, 8, 7, 0) in seen}")

1 + ... + n = n(n+1)/2 for n = 1..1000: True;  sum of squares = n(n+1)(2n+1)/6: True
factor(360) = [2, 2, 2, 3, 3, 5] in 11 calls; factor(9973) = [9973]
postage: 16c = 4x4 + 0x5, 23c = 2x4 + 3x5, 31c = 4x4 + 3x5, 100c = 25x4 + 0x5
3^1000 by repeated squaring: 16 multiplications instead of 999, invariant held at every step, correct: True
8-puzzle: 181,440 of 362,880 arrangements reachable from the goal; all with an even number of inversions: True
swapping tiles 7 and 8 gives 1 inversion: reachable: False

Invariants

An invariant is a property that holds initially and is preserved by every step. By induction on the number of steps, it then holds after any number of steps. This is the main tool for proving programs correct, an idea due to Robert Floyd and Tony Hoare in the late 1960s. For a loop, you show three things: the invariant holds before the first iteration; if it holds before an iteration, it holds after it; and when the loop exits, the invariant plus the exit condition imply the result you wanted.

Fast exponentiation keeps the invariant $r\cdot b^e=x^n$, where r is the result so far, b the base and e the remaining exponent. Initially r = 1, b = x and e = n. Each iteration either multiplies r by b when e is odd, or squares b and halves e, and either way the product is unchanged. When e reaches 0, the invariant says r = xⁿ. The example asserts the invariant on every iteration and computes 3¹⁰⁰⁰ in 16 multiplications instead of 999, because e halves each time. That's chapter 6's modular exponentiation, the operation at the heart of RSA.

Invariants also prove that something is impossible. In the 8-puzzle, read the eight tiles row by row and count the inversions, pairs that are out of order. Sliding a tile sideways doesn't change the reading order at all. Sliding a tile up or down moves it past exactly two other tiles, which changes the inversion count by −2, 0 or +2. Either way the parity of the count never changes. The solved position has 0 inversions, so no position with an odd number can ever reach it. Swapping tiles 7 and 8 creates one inversion, and no sequence of moves will ever undo it. The example's search confirms both halves of the story. Exactly 181,440 positions, half of the 9! arrangements, are reachable, and all of them have even parity, so parity is the only obstacle. The 15-puzzle craze of 1880 made a famous challenge of undoing a swap of 14 and 15 on the 4 × 4 board, and an invariant of the same kind shows that nobody ever could.

When induction goes wrong

The classic fake proof claims that all horses are the same colour. Base case: one horse is the same colour as itself. Step: given n + 1 horses, the first n are the same colour by hypothesis, and so are the last n, and since the two groups overlap, all n + 1 are. The step is fine for n ≥ 2, but from n = 1 to n = 2 the "first one" and the "last one" don't overlap at all, so the chain breaks at its first link. Check that the inductive step works for every n, including the smallest, and especially when an argument talks about "the others" or "the overlap".

Sets

A set is an unordered collection of distinct elements. Write $x\in A$ for membership, $A\subseteq B$ when every element of A is in B, and use union $A\cup B$, intersection $A\cap B$ and difference $A\setminus B$. The power set $𝒫(A)$ is the set of all subsets of A, and a set of n elements has $2^n$ of them, since each element is either in or out. The Cartesian product $A\times B$ is the set of ordered pairs (a, b). Set-builder notation, $\{x\in A\mid P(x)\}$, is a set comprehension: Python's {x for x in A if p(x)} is the same idea with the same shape.

Relations and their properties

A relation R on a set A is a set of ordered pairs from $A\times A$. Write a R b when (a, b) is one of them, like a ≤ b or "a divides b". Four properties classify relations:

• reflexive: a R a for every a;
• symmetric: a R b implies b R a;
• antisymmetric: a R b and b R a imply a = b;
• transitive: a R b and b R c imply a R c.

Each property is a "for all" statement over a finite set, so on a small set a program can check it exhaustively, as chapter 1 did with equivalences. Symmetric and antisymmetric aren't opposites. Equality is both, and "differs by at most 1" is symmetric without being transitive: 1 is close to 2 and 2 is close to 3, but 1 isn't close to 3. That relation is the trap behind any "approximately equal" comparison used as if it were equality.

Equivalence relations

A relation that is reflexive, symmetric and transitive is an equivalence relation, a precise version of "the same in the respect that matters". Every equivalence relation splits its set into disjoint equivalence classes, and every split into disjoint parts comes from an equivalence relation. Congruence mod 3 splits the integers into three classes by remainder, which is the foundation of modular arithmetic in chapter 5. "In the same connected component" is an equivalence relation on the vertices of a graph, and union-find is a data structure for maintaining its classes as edges arrive.

Programming languages quietly require equivalence relations. A hash table assumes that key equality is reflexive, symmetric and transitive, and that equal keys have equal hashes. Floating-point NaN breaks reflexivity, because nan == nan is false. In Python, two separately created NaN keys therefore land in a dictionary as two different keys, and you can't look up either one with a fresh NaN. A custom __eq__ that compares within a tolerance breaks transitivity, with similarly strange results.

Partial orders

A relation that is reflexive, antisymmetric and transitive is a partial order. It's "partial" because some pairs may be incomparable: neither 4 divides 6 nor 6 divides 4. Dependencies form partial orders, such as tasks in a build, courses with prerequisites and package versions, and so does "happened before" between events in a distributed system. A partial order is drawn as a Hasse diagram, with smaller elements lower and only the direct relationships shown, so that the rest follow by transitivity.

A chain is a set of mutually comparable elements, such as 1, 2, 4, 8. An antichain is a set of mutually incomparable ones, such as 7 to 12, none of which divides another. These two ideas answer scheduling questions. If the order is a set of dependencies, tasks in an antichain can run in parallel, and the longest chain is the critical path. You need at least as many rounds as the longest chain has elements, and grouping tasks by the length of the longest chain below them achieves exactly that number (Mirsky's theorem). Every finite partial order can also be listed in an order consistent with it, a topological sort, which is how build tools and package managers decide what to do first (exercise 4.2).

Functions

A function $f:A\to B$ assigns exactly one element of B to each element of A. It is injective (one-to-one) if different inputs give different outputs, surjective (onto) if every element of B is hit, and bijective if both. From an n-set to an m-set there are $m^n$ functions, because each of the n inputs independently picks one of m outputs. There are $m(m-1)\cdots (m-n+1)$ injective ones, because each input must avoid the outputs already used. A bijection between two finite sets proves they have the same size without counting either. That idea, the bijection rule, drives much of chapter 9, and the fact that no injection goes from a larger set to a smaller one is the pigeonhole principle of chapter 10.

Checked by enumeration

The example classifies four relations on 1 to 12, lists the classes of congruence mod 3, finds a longest chain and an antichain under divisibility, counts functions between small sets by brute force, and runs Cantor's diagonal construction on a finite list.

import random
from itertools import product

A = range(1, 13)
relations = {
    "a divides b":         lambda a, b: b % a == 0,
    "a = b (mod 3)":       lambda a, b: (a - b) % 3 == 0,
    "a <= b":              lambda a, b: a <= b,
    "|a - b| <= 1":        lambda a, b: abs(a - b) <= 1,
}
def props(R):
    refl = all(R(a, a) for a in A)
    sym = all(R(b, a) for a in A for b in A if R(a, b))
    anti = all(a == b for a in A for b in A if R(a, b) and R(b, a))
    trans = all(R(a, c) for a in A for b in A for c in A if R(a, b) and R(b, c))
    kind = ("equivalence relation" if refl and sym and trans else
            "partial order" if refl and anti and trans else "neither")
    return refl, sym, anti, trans, kind
print(f"{'relation on 1..12':18} reflexive symmetric antisymmetric transitive")
for name, R in relations.items():
    r, s, a, t, kind = props(R)
    print(f"{name:18} {r!s:9} {s!s:9} {a!s:13} {t!s:10} -> {kind}")

classes = {}
for a in A: classes.setdefault(a % 3, []).append(a)
print("equivalence classes of mod 3:", list(classes.values()))

# The divisibility order: chains (all comparable) and antichains (none comparable)
div = relations["a divides b"]
longest = {}
for b in A:                                    # longest chain ending at b, by dynamic programming
    longest[b] = max((longest[a] + [b] for a in A if a < b and div(a, b)), key=len, default=[b])
chain = max(longest.values(), key=len)
antichain = [a for a in A if a > 6]            # nothing in 7..12 divides anything else in 7..12
assert not any(div(a, b) for a in antichain for b in antichain if a != b)
print(f"longest chain under 'divides': {chain};  an antichain of size {len(antichain)}: {antichain}")

# Counting functions between finite sets by brute force
def count(n, m):
    fs = list(product(range(m), repeat=n))     # a function from n items to m items is a tuple
    inj = sum(len(set(f)) == n for f in fs); sur = sum(len(set(f)) == m for f in fs)
    return len(fs), inj, sur
for n, m in ((3, 4), (4, 3), (4, 4)):
    total, inj, sur = count(n, m)
    print(f"functions from a {n}-set to a {m}-set: {total:3} total, {inj:2} one-to-one, {sur:2} onto")

# Cantor's diagonal argument, run on a finite list
random.seed(4)
rows = ["".join(random.choice("01") for _ in range(8)) for _ in range(8)]
diag = "".join("1" if rows[i][i] == "0" else "0" for i in range(8))
print("listed:", " ".join(rows))
print(f"diagonal flipped: {diag}  in the list: {diag in rows} (it differs from row i at position i)")

relation on 1..12  reflexive symmetric antisymmetric transitive
a divides b        True      False     True          True       -> partial order
a = b (mod 3)      True      True      False         True       -> equivalence relation
a <= b             True      False     True          True       -> partial order
|a - b| <= 1       True      True      False         False      -> neither
equivalence classes of mod 3: [[1, 4, 7, 10], [2, 5, 8, 11], [3, 6, 9, 12]]
longest chain under 'divides': [1, 2, 4, 8];  an antichain of size 6: [7, 8, 9, 10, 11, 12]
functions from a 3-set to a 4-set:  64 total, 24 one-to-one,  0 onto
functions from a 4-set to a 3-set:  81 total,  0 one-to-one, 36 onto
functions from a 4-set to a 4-set: 256 total, 24 one-to-one, 24 onto
listed: 01011000 01100110 01001100 11101100 01101011 01111100 11000111 10001011
diagonal flipped: 10110000  in the list: False (it differs from row i at position i)

The table matches the definitions: divisibility and ≤ are partial orders, congruence mod 3 is an equivalence relation, and "within 1" is neither, since it fails transitivity. The function counts agree with the formulas: $4^3=64$ functions from a 3-set to a 4-set, of which 4 × 3 × 2 = 24 are one-to-one and none are onto, because three inputs can't cover four outputs. From a 4-set to a 4-set, the one-to-one and onto functions are exactly the same 24, the 4! permutations. For finite sets of equal size, injective, surjective and bijective coincide.

Countable and uncountable

Two sets have the same cardinality if there's a bijection between them. For infinite sets that has surprising consequences. The even numbers are as numerous as all natural numbers, via n ↔ 2n. The integers are too, listed 0, 1, −1, 2, −2, and so on, and even the rationals can be listed in a zig-zag through the grid of fractions p/q. A set that can be listed this way is countable. The set of all finite strings over a finite alphabet is countable: list them by length, then alphabetically. Every program is a finite string, so there are countably many programs.

Georg Cantor showed in 1891 that the infinite binary sequences are not countable. Given any list of them, build a new sequence whose i-th bit is the opposite of the i-th bit of the i-th sequence. The new sequence differs from every listed sequence in at least one position, so the list was incomplete. The example runs this on eight rows of eight bits, and the flipped diagonal 10110000 is indeed missing. Now read an infinite binary sequence as a yes-or-no problem about natural numbers, where bit n answers "is n in the set?". There are uncountably many such problems and only countably many programs, so almost every problem has no program that solves it. Turing's proof that the halting problem is undecidable uses the same diagonal trick to name one specific such problem.

Divisibility and primes

An integer a divides b, written $a\mid b$, if b = ak for some integer k. For any integer a and positive b, the division theorem gives a unique quotient q and remainder r with

and r is what a % b computes. A prime is an integer p > 1 whose only positive divisors are 1 and p. Chapter 3 proved that every integer above 1 is a product of primes. The fundamental theorem of arithmetic adds that the product is unique apart from order, and the proof of uniqueness needs a fact about primes proved below, Euclid's lemma.

Primes thin out, but slowly. The prime number theorem, proved in 1896, says that the number of primes up to n, written $\pi (n)$, satisfies $\pi (n)\sim n/\ln n$: the ratio tends to 1. So near n, roughly one number in $\ln n$ is prime. Near $2^{1024}$ that's one in about 710, or one odd number in about 355, which is why generating a random 1024-bit prime takes a few hundred tries rather than billions. The example computes $\pi (n)$ up to ten million with a sieve, and the ratio creeps towards 1: 1.161 at a thousand, 1.071 at ten million.

Euclid's algorithm

The greatest common divisor gcd(a, b) is the largest integer dividing both. Euclid's algorithm, from around 300 BC, rests on one identity:

Proof. Write a = qb + r. Any common divisor of a and b divides r = a − qb, and any common divisor of b and r divides a = qb + r. So the two pairs have exactly the same common divisors, and in particular the same greatest one. ∎ Repeating the identity shrinks the second number each time, and when it reaches 0, gcd(g, 0) = g is the answer. Geometrically, each step cuts the largest possible squares off a rectangle, and the side of the last square is the gcd.

How many steps can it take? The remainders shrink fastest when quotients are large and slowest when every quotient is 1, which happens exactly for consecutive Fibonacci numbers. Gabriel Lamé proved in 1844 that the number of steps is at most about five times the number of decimal digits of the smaller number, so Euclid's algorithm is logarithmic. This was one of the first results in the analysis of algorithms. The example runs the worst case, two Fibonacci numbers near a million that take 29 steps, while random pairs of the same size average under 12.

Bézout's identity and modular inverses

Running Euclid's algorithm backwards gives more than the gcd. Bézout's identity says there are integers x and y with

and the extended Euclidean algorithm finds them by carrying the coefficients through each step. For example, gcd(240, 46) = 2 = 240 × (−9) + 46 × 47. Two consequences follow at once.

Inverses. If gcd(a, n) = 1, then ax + ny = 1, so $ax\equiv 1(\mathrm{mod}n)$: x is an inverse of a modulo n. If gcd(a, n) = d > 1, then every multiple of a mod n is a multiple of d, so none equals 1 and no inverse exists. The example lists the inverses mod 7, where every nonzero number has one, and mod 8, where only the odd numbers do. Python computes inverses directly with pow(a, -1, n).

Euclid's lemma. If a prime p divides ab, then p divides a or p divides b. Proof. Suppose $p\nmid a$. Then gcd(a, p) = 1, so 1 = ax + py for some x, y. Multiply by b: b = abx + pby. Both terms are multiples of p, so $p\mid b$. ∎ Uniqueness of prime factorisation follows: if two factorisations differed, a prime from one would divide the product of the other, and so, by the lemma, one of its primes, which then must be the same prime.

Modular arithmetic

Chapter 4 showed that congruence mod n, $a\equiv b(\mathrm{mod}n)$, is an equivalence relation. It also respects arithmetic: if $a\equiv a^{\prime }$ and $b\equiv b^{\prime }$, then $a+b\equiv a^{\prime }+b^{\prime }$ and $ab\equiv a^{\prime }{b}^{\prime }$. So you can reduce mod n at any point in a calculation and the final remainder doesn't change. That's what keeps the numbers in chapter 6 small enough to handle, and it's how hash functions, checksums and clock arithmetic work. Addition, subtraction and multiplication behave as usual. Division is the one difference: you can divide by a only if a has an inverse mod n. When n is a prime p, every nonzero residue has one, so arithmetic mod p supports all four operations, a structure called a finite field.

One practical trap: Python's % always returns a result with the sign of the divisor, so -7 % 3 is 2, but C, Java, Go and JavaScript return −1. Code that indexes a hash table with h % size can produce negative indexes in those languages when h is negative.

The Chinese remainder theorem

A problem in the Sunzi Suanjing, a Chinese mathematical text from around the 3rd to 5th century, asks for a number that leaves remainder 2 when divided by 3, remainder 3 when divided by 5, and remainder 2 when divided by 7. The Chinese remainder theorem says that if the moduli $m_1,\dots ,m_k$ are pairwise coprime, such a system has exactly one solution modulo their product $M=m_1{m}_2\cdots m_k$. The proof is a construction. With $M_i=M/m_i$, the number

works, because the i-th term is $r_i$ mod $m_i$ and 0 mod every other modulus. The inverse exists because $M_i$ is coprime to $m_i$. Uniqueness holds because two solutions differ by a multiple of every $m_i$, hence of M. The answer to the old problem is 23. In computing, the theorem lets you replace one computation mod a large M with several mod small factors. RSA implementations decrypt mod p and mod q separately and recombine, which is about four times faster than working mod n directly.

All of it, computed

The example sieves the primes up to ten million, runs Euclid's algorithm on its worst case and on random pairs, finds Bézout coefficients and inverses, and solves the Sunzi problem by the construction above.

import math, random

# 1. The sieve of Eratosthenes, and how many primes there are
N = 10 ** 7
sieve = bytearray([1]) * (N + 1); sieve[0] = sieve[1] = 0
for p in range(2, math.isqrt(N) + 1):
    if sieve[p]: sieve[p * p::p] = bytearray(len(range(p * p, N + 1, p)))
print("     n      primes up to n   n / ln n   ratio")
count = 0; marks = {10 ** k for k in range(1, 8)}
for n in range(N + 1):
    count += sieve[n]
    if n in marks: print(f"{n:>10,}  {count:>12,}  {n / math.log(n):>11,.0f}   {count / (n / math.log(n)):.3f}")

# 2. Euclid's algorithm, counting division steps
def gcd_steps(a, b):
    steps = 0
    while b: a, b = b, a % b; steps += 1
    return a, steps
fib = [1, 1]
while len(fib) < 31: fib.append(fib[-1] + fib[-2])
print(f"gcd({fib[30]:,}, {fib[29]:,}) = {gcd_steps(fib[30], fib[29])[0]} after {gcd_steps(fib[30], fib[29])[1]} steps "
      "(consecutive Fibonacci numbers are the worst case)")
random.seed(5)
pairs = [(random.randrange(1, 10 ** 6), random.randrange(1, 10 ** 6)) for _ in range(100_000)]
avg = sum(gcd_steps(a, b)[1] for a, b in pairs) / len(pairs)
print(f"average steps for 100,000 random pairs below 1,000,000: {avg:.2f}")

# 3. Extended Euclid: Bezout coefficients and modular inverses
def ext_gcd(a, b):
    if b == 0: return a, 1, 0
    g, x, y = ext_gcd(b, a % b)
    return g, y, x - (a // b) * y              # from b*x + (a mod b)*y = g
g, x, y = ext_gcd(240, 46)
print(f"gcd(240, 46) = {g} = 240 x ({x}) + 46 x {y}")
g, x, _ = ext_gcd(17, 3120)
print(f"inverse of 17 mod 3120: {x % 3120}  (check: 17 x {x % 3120} mod 3120 = {17 * (x % 3120) % 3120});  "
      f"Python: pow(17, -1, 3120) = {pow(17, -1, 3120)}")
for m in (7, 8):
    inv = {a: next((b for b in range(1, m) if a * b % m == 1), None) for a in range(1, m)}
    print(f"inverses mod {m}: " + ", ".join(f"{a}->{b if b else '-'}" for a, b in inv.items()))

# 4. The Chinese remainder theorem: x = 2 mod 3, x = 3 mod 5, x = 2 mod 7
def crt(residues, moduli):
    M = math.prod(moduli); x = 0
    for r, m in zip(residues, moduli):
        Mi = M // m
        x += r * Mi * pow(Mi, -1, m)           # this term is r mod m and 0 mod every other modulus
    return x % M, M
x, M = crt([2, 3, 2], [3, 5, 7])
print(f"x = 2 (mod 3), 3 (mod 5), 2 (mod 7)  =>  x = {x} (mod {M});  check: {x % 3}, {x % 5}, {x % 7}")

     n      primes up to n   n / ln n   ratio
        10             4            4   0.921
       100            25           22   1.151
     1,000           168          145   1.161
    10,000         1,229        1,086   1.132
   100,000         9,592        8,686   1.104
 1,000,000        78,498       72,382   1.084
10,000,000       664,579      620,421   1.071
gcd(1,346,269, 832,040) = 1 after 29 steps (consecutive Fibonacci numbers are the worst case)
average steps for 100,000 random pairs below 1,000,000: 11.69
gcd(240, 46) = 2 = 240 x (-9) + 46 x 47
inverse of 17 mod 3120: 2753  (check: 17 x 2753 mod 3120 = 1);  Python: pow(17, -1, 3120) = 2753
inverses mod 7: 1->1, 2->4, 3->5, 4->2, 5->3, 6->6
inverses mod 8: 1->1, 2->-, 3->3, 4->-, 5->5, 6->-, 7->7
x = 2 (mod 3), 3 (mod 5), 2 (mod 7)  =>  x = 23 (mod 105);  check: 2, 3, 2

The inverse of 17 mod 3120 is 2753, which is about to become a private key: chapter 6's tiny RSA example uses exactly these numbers. The inverse tables show the rule about division in action. Mod 7, every nonzero number has an inverse. Mod 8, the even numbers have none, because they share the factor 2 with the modulus.

Fermat's little theorem

For a prime p and any integer a not divisible by p,

Proof. Multiply each of the numbers 1, 2, …, p − 1 by a, mod p. The results are all nonzero, since p divides neither factor, and all different, since a has an inverse mod p and $ai\equiv aj$ would give $i\equiv j$. So the products are just 1, …, p − 1 in some order, and multiplying them all together gives

Every factor of (p − 1)! is invertible mod p, so cancel it to get $a^{p-1}\equiv 1$. ∎ Leonhard Euler generalised the theorem to any modulus. Let $\phi (n)$ count the numbers from 1 to n that are coprime to n. Then $a^{\phi (n)}\equiv 1(\mathrm{mod}n)$ whenever gcd(a, n) = 1, by the same argument applied to those $\phi (n)$ numbers. For n = pq with distinct primes, $\phi (n)=(p-1)(q-1)$, since the numbers not coprime to n are the q multiples of p and the p multiples of q, which share n itself.

RSA

Ron Rivest, Adi Shamir and Leonard Adleman published the scheme in 1977. To make a key pair:

• pick two large random primes p and q, and let n = pq;
• pick a public exponent e coprime to (p − 1)(q − 1), almost always 65537;
• compute the private exponent $d=e^{-1}\mathrm{mod}\phi (n)$ with the extended Euclidean algorithm. Implementations often use $\lambda (n)=\mathrm{lcm}(p-1,q-1)$ instead of $\phi (n)$, which gives a smaller d that works the same way.

The public key is (n, e) and the private key is d. To encrypt a message m < n, compute $c=m^e\mathrm{mod}n$. To decrypt, compute $c^d\mathrm{mod}n$. Both are fast thanks to chapter 3's square-and-multiply with a reduction mod n after each step.

Why decryption works. Since $ed\equiv 1(\mathrm{mod}\phi (n))$, write $ed=1+k\phi (n)$. If gcd(m, n) = 1, Euler's theorem gives

The rare case where m shares a factor with n needs one more step through the Chinese remainder theorem (exercise 6.2). The same algebra gives signatures: Alice computes $s=m^d\mathrm{mod}n$, and anyone can check $s^e\mathrm{mod}n=m$. Only the holder of d could have produced s. With chapter 5's numbers, p = 61 and q = 53 give n = 3233, e = 17 and d = 2753, and the message 65 encrypts to 2790 and decrypts back to 65.

Finding large primes

Key generation needs random primes of 1024 bits or more. By the prime number theorem, about one odd number in 355 of that size is prime, so the generator draws random odd candidates and tests each. It needs a fast test, because trial division up to $\sqrt{n}$ would take longer than the age of the universe.

Fermat's little theorem suggests one: if $a^{n-1}\not\equiv 1(\mathrm{mod}n)$ for some a, then n is certainly composite. The converse fails, though. Some composites, the Carmichael numbers, satisfy $a^{n-1}\equiv 1$ for every a coprime to n. The smallest is 561 = 3 × 11 × 17, and in 1994 Alford, Granville and Pomerance proved there are infinitely many. The Miller–Rabin test closes the gap. Write $n-1=2^{s}d$ with d odd. For a random base a, look at the sequence $a^d,a^{2d},\dots ,a^{2^{s}d}$ mod n. If n is prime, the only square roots of 1 mod n are $\pm 1$, so the sequence must either start at 1 or pass through $-1$ before reaching 1. A composite n violates this for at least three quarters of all bases, so 40 independent rounds leave an error probability below $4^{-40}\approx {10}^{-24}$. That's far less likely than a hardware fault during the computation.

Built from scratch

The example finds every Carmichael number below 10,000 and shows Miller–Rabin rejecting each one, repeats the tiny RSA example, generates two random 1024-bit primes, and uses the 2048-bit key to encrypt, decrypt, sign and verify. It then demonstrates two weaknesses of textbook RSA. Its random numbers come from Python's random with a fixed seed so that the output is reproducible, which makes these keys worthless. Real keys must come from a cryptographically secure source such as Python's secrets module or the operating system.

import math, random

def is_probable_prime(n, rounds=40, rng=random):
    """Miller-Rabin: write n - 1 = 2^s * d; a witness a exposes a composite n."""
    if n < 4: return n in (2, 3)
    if n % 2 == 0: return False
    s, d = 0, n - 1
    while d % 2 == 0: s += 1; d //= 2
    for _ in range(rounds):
        a = rng.randrange(2, n - 1)
        x = pow(a, d, n)
        if x in (1, n - 1): continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1: break
        else:
            return False                       # a is a witness: n is certainly composite
    return True                                # no witness found: n is prime with high probability

# 1. Fermat's test is fooled by Carmichael numbers; Miller-Rabin is not
def fermat_passes(n): return all(pow(a, n - 1, n) == 1 for a in range(2, n) if math.gcd(a, n) == 1)
def composite(n): return any(n % d == 0 for d in range(2, math.isqrt(n) + 1))
carmichael = [n for n in range(3, 10_000, 2) if composite(n) and fermat_passes(n)]
random.seed(6)
print(f"composites below 10,000 that pass Fermat's test for every coprime base: {carmichael}")
print(f"561 = 3 x 11 x 17; Miller-Rabin verdicts: {[is_probable_prime(n) for n in carmichael]}")

# 2. Textbook RSA with tiny numbers
p, q, e = 61, 53, 17
n, phi = p * q, (p - 1) * (q - 1)
d = pow(e, -1, phi)
c = pow(65, e, n)
print(f"\ntiny RSA: n = {n}, e = {e}, d = {d};  65 -> {c} -> {pow(c, d, n)}")

# 3. Real-sized keys: random 1024-bit primes
def random_prime(bits):
    tries = 0
    while True:
        tries += 1
        cand = random.getrandbits(bits) | (1 << (bits - 1)) | 1     # top bit set, odd
        if is_probable_prime(cand): return cand, tries
p, tp = random_prime(1024); q, tq = random_prime(1024)
n, e = p * q, 65537
lam = math.lcm(p - 1, q - 1)
d = pow(e, -1, lam)
print(f"1024-bit primes found after {tp} and {tq} odd candidates (expected about {1024 * math.log(2) / 2:.0f} each)")
m = int.from_bytes(b"attack at dawn", "big")
c = pow(m, e, n)
back = pow(c, d, n).to_bytes(14, "big")
print(f"n has {n.bit_length()} bits; decrypt(encrypt('attack at dawn')) = {back.decode()!r}")
sig = pow(m, d, n)
print(f"signature verifies: {pow(sig, e, n) == m};  altered message verifies: {pow(sig, e, n) == m + 1}")

# 4. Why textbook RSA is unsafe without padding
c2 = c * pow(2, e, n) % n                      # anyone can do this without the private key
print(f"multiplying the ciphertext by 2^e mod n doubles the plaintext: {pow(c2, d, n) == 2 * m}")
print(f"encrypting the same message twice gives the same ciphertext: {pow(m, e, n) == c}")

# 5. Cost: square-and-multiply on a 2048-bit exponent
squarings, multiplies = d.bit_length() - 1, bin(d).count("1") - 1
print(f"decryption exponent has {d.bit_length()} bits: about {squarings} squarings + {multiplies} multiplications")

composites below 10,000 that pass Fermat's test for every coprime base: [561, 1105, 1729, 2465, 2821, 6601, 8911]
561 = 3 x 11 x 17; Miller-Rabin verdicts: [False, False, False, False, False, False, False]

tiny RSA: n = 3233, e = 17, d = 2753;  65 -> 2790 -> 65
1024-bit primes found after 55 and 275 odd candidates (expected about 355 each)
n has 2048 bits; decrypt(encrypt('attack at dawn')) = 'attack at dawn'
signature verifies: True;  altered message verifies: False
multiplying the ciphertext by 2^e mod n doubles the plaintext: True
encrypting the same message twice gives the same ciphertext: True
decryption exponent has 2047 bits: about 2046 squarings + 995 multiplications

The two primes took 55 and 275 odd candidates, around the expected 355, because the number of tries before a success is geometrically distributed (chapter 14). Decryption with a 2047-bit exponent costs about 2,046 squarings and 995 multiplications of 2048-bit numbers, which is why RSA private-key operations are expensive and why implementations use the Chinese remainder theorem to speed them up.

Why textbook RSA is broken

Everything above is mathematically correct, and the last two lines of output show that it's still insecure. First, it's deterministic. The same message always gives the same ciphertext, so an attacker who can guess the possible messages, "yes" or "no", a salary from a short list, simply encrypts each guess with the public key and compares. Second, it's malleable. Multiplying a ciphertext by $2^e$ doubles the hidden plaintext without anyone knowing it, and real attacks have exploited structure like this. Small exponents with small messages are worse still: with e = 3 and $m^3<n$, the ciphertext is just $m^3$ and an ordinary cube root recovers m. Real RSA therefore always uses randomised padding, OAEP for encryption and PSS for signatures, and modern protocols mostly use elliptic-curve Diffie–Hellman for key exchange. The Engineer's Citadel covers what to use in practice. The rule from this chapter is simpler: understand the mathematics, then never implement it yourself in production.

The security also rests on an assumption. Nobody has proved factoring is hard. The largest RSA challenge number factored so far, RSA-250 with 829 bits, fell in 2020 after thousands of core-years of computation, and 2048-bit keys are the current minimum recommendation. A large quantum computer running Shor's algorithm would factor efficiently, which is why NIST standardised the post-quantum key-encapsulation scheme ML-KEM in 2024.

Vocabulary

A graph G = (V, E) is a set of vertices V and a set of edges E, each edge joining two vertices. In a simple graph there are no loops and no repeated edges; a multigraph allows repeats. The degree of a vertex is the number of edges at it. A walk is a sequence of vertices joined by edges, a path is a walk that repeats no vertex, and a cycle is a closed walk that repeats no vertex except the start. A graph is connected if there's a path between every pair of vertices, and its maximal connected pieces are its components, the classes of the equivalence relation "is connected to" from chapter 4. Directed graphs give each edge a direction, and a vertex then has separate in- and out-degrees.

The handshake lemma

In any graph,

Proof. Each edge has two ends and contributes exactly 1 to the degree of each, so summing degrees counts every edge twice. ∎ This is double counting, counting the same thing in two ways, a technique chapter 9 uses again. One corollary follows at once. The sum of all degrees is even, so the number of vertices with odd degree must be even: in any group, the number of people who have shaken an odd number of hands is even. The example's random graph has 993 edges, degree sum 1986, and 96 odd-degree vertices, and that last number will be even whatever the random seed.

Euler tours

Graph theory began with a puzzle about the city of Königsberg, whose seven bridges joined two river banks and two islands. Could you walk through the city crossing every bridge exactly once? Leonhard Euler showed in 1736 that you couldn't. Model each land mass as a vertex and each bridge as an edge.

Euler's theorem. A connected graph has a closed walk that uses every edge exactly once, an Euler circuit, if and only if every vertex has even degree. It has an open walk that does so if and only if exactly two vertices have odd degree, and then the walk must start at one and end at the other. Necessity is the argument in the caption. Sufficiency is constructive. Start anywhere and walk along unused edges. Because degrees are even, you can only get stuck back at the start. If unused edges remain, some vertex on your circuit touches one, so walk a detour from there and splice it in. That's Hierholzer's algorithm, linear in the number of edges, and the example uses it to find a circuit through all ten edges of K₅, where every degree is 4.

Change one word and everything changes. A Hamiltonian cycle visits every vertex exactly once, and deciding whether one exists is NP-complete (The Engineer's Workbench, chapter 8). No simple degree condition characterises it. The contrast is a good example of how a problem's statement, not its look, decides whether it's easy.

Trees

A tree is a connected graph with no cycles. For a graph with n vertices, the following are equivalent, and each is a useful definition:

• it's connected and has no cycles;
• it's connected and has exactly n − 1 edges;
• it has no cycles and has exactly n − 1 edges;
• between any two vertices there's exactly one path.

The key step is that every tree with n ≥ 2 vertices has a leaf, a vertex of degree 1: take a longest path, and its end can't have any other neighbour, since that would extend the path or close a cycle. Removing a leaf and its edge leaves a tree on n − 1 vertices, so by induction a tree on n vertices has n − 1 edges. The example's spanning-tree test relies on the third form: n − 1 edges that never close a cycle must form a tree. It counts the spanning trees of K_n by brute force and matches Cayley's formula, $n^{n-2}$: 16, 125 and 1296. Spanning trees appear in network protocols that switch off redundant links to avoid loops, and in the minimum spanning tree algorithms of stage 02.

Bipartite graphs

A graph is bipartite if its vertices split into two sides with every edge crossing between them: workers and jobs, students and courses, users and items. Theorem: a graph is bipartite if and only if it has no cycle of odd length. One direction is easy. A cycle alternates sides with each step, so it returns to its start only after an even number of steps. For the other, run breadth-first search from any vertex and colour vertices by the parity of their distance. If some edge joins two vertices of the same colour, the two BFS paths to its ends plus the edge itself form a closed walk of odd length, which contains an odd cycle. So if there's no odd cycle, the colouring is a valid split. The proof is the algorithm: BFS either 2-colours the graph or exposes an odd cycle. In the example, a 6-cycle and an 8 × 8 grid are bipartite (colour the grid like a chessboard), while a 7-cycle and the Petersen graph, which contains 5-cycles, are not.

Checked, constructed and simulated

The example verifies the handshake lemma, applies Euler's degree test to Königsberg and constructs a circuit in K₅, counts spanning trees, tests four graphs for bipartiteness, and then builds random graphs on 20,000 vertices with increasing average degree.

import math, random
from collections import deque, Counter
from itertools import combinations

random.seed(7)
# 1. The handshake lemma on a random graph
n = 200
edges = [(a, b) for a, b in combinations(range(n), 2) if random.random() < 0.05]
deg = Counter(v for e in edges for v in e)
odd = sum(deg[v] % 2 for v in range(n))
print(f"random graph: {n} vertices, {len(edges)} edges; sum of degrees = {sum(deg.values())} = 2 x {len(edges)}; "
      f"{odd} vertices of odd degree (always even)")

# 2. Euler: the bridges of Konigsberg, and a circuit through every edge of K5
bridges = [("A", "B"), ("A", "B"), ("A", "C"), ("A", "C"), ("A", "D"), ("B", "D"), ("C", "D")]
kdeg = Counter(v for e in bridges for v in e)
print(f"Konigsberg degrees {dict(sorted(kdeg.items()))}: {sum(d % 2 for d in kdeg.values())} odd vertices, "
      "so no walk crosses every bridge exactly once")
def euler_circuit(edge_list, start):
    adj = {}
    for i, (a, b) in enumerate(edge_list):
        adj.setdefault(a, []).append((b, i)); adj.setdefault(b, []).append((a, i))
    used, stack, circuit = set(), [start], []
    while stack:                                # Hierholzer: walk until stuck, then splice in detours
        v = stack[-1]
        while adj[v] and adj[v][-1][1] in used: adj[v].pop()
        if adj[v]:
            w, i = adj[v].pop(); used.add(i); stack.append(w)
        else:
            circuit.append(stack.pop())
    return circuit[::-1]
k5 = list(combinations(range(5), 2))
tour = euler_circuit(k5, 0)
steps = {frozenset(p) for p in zip(tour, tour[1:])}
print(f"K5 (every degree 4): Euler circuit {'-'.join(map(str, tour))} uses all {len(steps)} edges once")

# 3. Trees: count the spanning trees of K_n by brute force, against Cayley's formula n^(n-2)
def is_tree(n, es):
    parent = list(range(n))
    def find(x):
        while parent[x] != x: parent[x] = parent[parent[x]]; x = parent[x]
        return x
    for a, b in es:
        ra, rb = find(a), find(b)
        if ra == rb: return False               # this edge would close a cycle
        parent[ra] = rb
    return True                                 # n - 1 edges and no cycle: connected, a tree
for m in (4, 5, 6):
    count = sum(is_tree(m, es) for es in combinations(list(combinations(range(m), 2)), m - 1))
    print(f"spanning trees of K{m}: {count:5,} by brute force, {m ** (m - 2):5,} by Cayley's formula")

# 4. Bipartite or not: BFS 2-colouring, which either succeeds or finds an odd cycle
def two_colour(adj):
    colour = {}
    for s in adj:
        if s in colour: continue
        colour[s] = 0; q = deque([s])
        while q:
            v = q.popleft()
            for w in adj[v]:
                if w not in colour: colour[w] = 1 - colour[v]; q.append(w)
                elif colour[w] == colour[v]: return False
    return True
def cycle(k): return {i: [(i - 1) % k, (i + 1) % k] for i in range(k)}
petersen = {i: [] for i in range(10)}
for i in range(5):
    for a, b in ((i, (i + 1) % 5), (i, i + 5), (i + 5, (i + 2) % 5 + 5)):
        petersen[a].append(b); petersen[b].append(a)
grid = {(r, c): [(r + dr, c + dc) for dr, dc in ((1, 0), (-1, 0), (0, 1), (0, -1)) if 0 <= r + dr < 8 and 0 <= c + dc < 8]
        for r in range(8) for c in range(8)}
print("2-colourable:", ", ".join(f"{name} {two_colour(g)}" for name, g in
      (("6-cycle", cycle(6)), ("7-cycle", cycle(7)), ("8x8 grid", grid), ("Petersen graph", petersen))))

# 5. Random graphs: a giant component appears suddenly when the average degree passes 1
N = 20_000
print(f"\naverage degree c, random graph on {N:,} vertices: share of vertices in the largest component")
for c in (0.5, 0.9, 1.0, 1.1, 1.5, 2.0, 3.0):
    parent = list(range(N))
    def find(x):
        while parent[x] != x: parent[x] = parent[parent[x]]; x = parent[x]
        return x
    for _ in range(round(c * N / 2)):
        a, b = random.randrange(N), random.randrange(N)
        parent[find(a)] = find(b)
    largest = max(Counter(find(v) for v in range(N)).values())
    beta = 0.0
    if c > 1:
        beta = 0.5
        for _ in range(200): beta = 1 - math.exp(-c * beta)   # fixed point of beta = 1 - e^(-c beta)
    print(f"   c = {c:3}: {largest / N:6.1%}   (theory for large n: {beta:.1%})")

random graph: 200 vertices, 993 edges; sum of degrees = 1986 = 2 x 993; 96 vertices of odd degree (always even)
Konigsberg degrees {'A': 5, 'B': 3, 'C': 3, 'D': 3}: 4 odd vertices, so no walk crosses every bridge exactly once
K5 (every degree 4): Euler circuit 0-4-3-2-4-1-3-0-2-1-0 uses all 10 edges once
spanning trees of K4:    16 by brute force,    16 by Cayley's formula
spanning trees of K5:   125 by brute force,   125 by Cayley's formula
spanning trees of K6: 1,296 by brute force, 1,296 by Cayley's formula
2-colourable: 6-cycle True, 7-cycle False, 8x8 grid True, Petersen graph False

average degree c, random graph on 20,000 vertices: share of vertices in the largest component
   c = 0.5:   0.1%   (theory for large n: 0.0%)
   c = 0.9:   1.2%   (theory for large n: 0.0%)
   c = 1.0:   3.9%   (theory for large n: 0.0%)
   c = 1.1:  18.1%   (theory for large n: 17.6%)
   c = 1.5:  58.3%   (theory for large n: 58.3%)
   c = 2.0:  79.4%   (theory for large n: 79.7%)
   c = 3.0:  94.1%   (theory for large n: 94.0%)

Random graphs and the giant component

The last table shows one of the most striking facts about networks. Paul Erdős and Alfréd Rényi studied graphs whose edges are chosen at random, and found a phase transition. When the average degree c is below 1, every component is small, a vanishing fraction of the graph: 0.1% at c = 0.5. Just above 1, a giant component emerges containing a fixed fraction $\beta$ of all vertices, where $\beta$ solves

At c = 1.5 the theory predicts 58.3% and the simulation gives 58.3%; at c = 3, 94.0% against 94.1%. The intuition is a branching process. Explore outward from a vertex: each newly reached vertex brings about c new neighbours on average. Below 1 the exploration dies out quickly. Above 1 it survives with positive probability, and the survivors all join one giant cluster. The same threshold appears in epidemics, where an outbreak becomes large when each case infects more than one other on average, and in the resilience of networks to random failures.

Colouring

A proper colouring gives each vertex a colour so that adjacent vertices differ, and the chromatic number $\chi (G)$ is the fewest colours that suffice. Colouring models any problem of conflicts. Exams that share a student need different slots, variables that are live at the same time need different registers, and nearby transmitters need different frequencies. Computing $\chi (G)$ is NP-hard, and even deciding whether three colours suffice is NP-complete. One bound is always easy, though.

Theorem: $\chi (G)\le \mathrm{\Delta }+1$, where $\mathrm{\Delta }$ is the maximum degree. Proof. Colour the vertices one at a time in any order, giving each the smallest colour not used by its already-coloured neighbours. A vertex has at most $\mathrm{\Delta }$ neighbours, so at most $\mathrm{\Delta }$ colours are blocked and one of the first $\mathrm{\Delta }+1$ is free. ∎ The bound says nothing about how good greedy is on a particular graph, and the order matters enormously. The example colours a crown graph, which is bipartite and so needs only 2 colours: in one order greedy uses 8, in another 2. On a random graph with maximum degree 29, greedy used 9 colours, far below the bound of 30, and ordering by degree first saved one more. That's typical. The bound is a guarantee, not a prediction.

Matching and Hall's theorem

A matching is a set of edges with no shared vertices. In a bipartite graph of people and jobs, it's an assignment where nobody gets two jobs and no job goes to two people. When can everyone on the left be matched? Write N(S) for the set of jobs that someone in a set S of people would accept. If some set S of people likes fewer than |S| jobs between them, they can't all be matched, by the pigeonhole principle. Hall's theorem (1935) says this is the only obstacle: a matching that covers the whole left side exists if and only if

The proof is an induction on the number of people, splitting on whether some proper subset S has |N(S)| exactly equal to |S|. Algorithms don't check all $2^n$ subsets. They grow a matching with augmenting paths, paths that alternate between unmatched and matched edges and start and end at unmatched vertices. Flipping such a path adds one edge to the matching, and when none exists the matching is maximum. When it falls short, the search also finds the set S that violates Hall's condition, a certificate of failure. In the example's four-person instance, the maximum matching covers only three people, and the reason is that ann, bob and cat between them accept only db and web, three people for two jobs.

Stable matching

Now suppose both sides have preferences: proposers rank receivers and receivers rank proposers. A matching is unstable if some proposer and receiver prefer each other to their assigned partners, because they'd both rather break away. David Gale and Lloyd Shapley showed in 1962 that a stable matching always exists, by an algorithm. Each free proposer proposes to the best receiver they haven't yet proposed to. Each receiver holds on to the best proposal so far and rejects the rest, trading up whenever a better one arrives.

It terminates: each proposer proposes to each receiver at most once, so there are at most $n^2$ proposals. It ends with everyone matched: a receiver, once proposed to, always holds someone, so a proposer rejected by everyone would mean all n receivers are held by the other n − 1 proposers, which is impossible. It's stable: suppose p prefers r to the partner p ended with. Then p proposed to r earlier and was rejected, or later dropped, in favour of someone r preferred, and since receivers only ever trade up, r prefers her final partner to p. So (p, r) doesn't block. ∎

The algorithm isn't neutral. Every proposer gets the best partner they could have in any stable matching, and every receiver gets the worst. The example measures it with 1,000 on each side and random preferences. Proposers end up with, on average, their 8th choice and receivers with their 123rd. Swap the roles with the same preferences, and the former proposers fall to an average of 155th. The National Resident Matching Program, which places US medical graduates in residencies, runs a variant of this algorithm, and in 2012 Shapley and Alvin Roth received the Nobel prize in economics for this work. After a long debate about who should propose, the program now runs it with applicants proposing.

Planarity

A graph is planar if it can be drawn in the plane with no edges crossing, which matters for laying out circuits on a single layer and drawing readable diagrams. A planar drawing divides the plane into faces, counting the unbounded outside, and Euler's formula connects the three counts. For any connected planar drawing,

Proof by induction on E. A connected graph with the fewest edges is a tree, with E = V − 1 and one face, and V − (V − 1) + 1 = 2. Adding an edge between two vertices of a connected drawing, without crossings, cuts one face into two: E and F each go up by 1, and the sum is unchanged. ∎ The example checks the cube, the octahedron, the icosahedron and a 10 × 10 grid.

The formula proves impossibility. In a simple planar graph with at least three vertices, every face is bounded by at least 3 edges and every edge borders at most 2 faces, so $3F\le 2E$. Substituting into Euler's formula gives

K₅ has 10 edges, but 3 × 5 − 6 = 9, so it can't be planar. If the graph has no triangles, every face needs at least 4 edges, which gives $E\le 2V-4$. K_3,3, the "three houses, three utilities" puzzle, has 9 edges against 2 × 6 − 4 = 8, so it's not planar either.

The same inequality gives planar graphs a vertex of degree at most 5. If every degree were 6 or more, the degree sum 2E would be at least 6V, contradicting $E\le 3V-6$. Remove that vertex, colour the rest by induction, and put it back. At most 5 neighbours block at most 5 colours, so every planar graph is 6-colourable. The famous four colour theorem brings this down to 4. Kenneth Appel and Wolfgang Haken proved it in 1976 with a computer checking nearly two thousand configurations, the first major theorem proved that way, and Georges Gonthier verified it entirely in the Coq proof assistant in 2005.

The algorithms, measured

The example compares greedy colouring orders, runs Gale–Shapley at scale and checks for blocking pairs, finds a maximum matching and the set that violates Hall's condition, and confirms Euler's formula and the edge bounds.

import random
from itertools import combinations

random.seed(8)
# 1. Greedy colouring: never more than (max degree + 1) colours, but the order matters
def greedy(order, adj):
    colour = {}
    for v in order:
        taken = {colour[w] for w in adj[v] if w in colour}
        colour[v] = next(c for c in range(len(adj) + 1) if c not in taken)
    return colour
k = 8                                           # crown graph: u_i joined to v_j whenever i != j
crown = {("u", i): [("v", j) for j in range(k) if j != i] for i in range(k)}
crown.update({("v", j): [("u", i) for i in range(k) if i != j] for j in range(k)})
bad = [x for i in range(k) for x in (("u", i), ("v", i))]
good = [("u", i) for i in range(k)] + [("v", i) for i in range(k)]
print(f"crown graph on {2 * k} vertices (bipartite): greedy uses {len(set(greedy(bad, crown).values()))} colours "
      f"in one order, {len(set(greedy(good, crown).values()))} in another")
n = 300
adj = {v: [] for v in range(n)}
for a, b in combinations(range(n), 2):
    if random.random() < 0.05: adj[a].append(b); adj[b].append(a)
maxdeg = max(len(x) for x in adj.values())
plain = len(set(greedy(range(n), adj).values()))
by_degree = len(set(greedy(sorted(adj, key=lambda v: -len(adj[v])), adj).values()))
print(f"random graph, max degree {maxdeg}: greedy uses {plain} colours (bound {maxdeg + 1}); largest-degree-first: {by_degree}")

# 2. Stable matching: Gale-Shapley with 1,000 proposers and 1,000 receivers
def gale_shapley(prop_pref, recv_pref):
    m = len(prop_pref)
    rank = [{p: r for r, p in enumerate(prefs)} for prefs in recv_pref]
    nxt, partner, free, proposals = [0] * m, [None] * m, list(range(m)), 0
    while free:
        p = free.pop()
        r = prop_pref[p][nxt[p]]; nxt[p] += 1; proposals += 1
        if partner[r] is None: partner[r] = p
        elif rank[r][p] < rank[r][partner[r]]: free.append(partner[r]); partner[r] = p
        else: free.append(p)
    return {partner[r]: r for r in range(m)}, proposals
m = 1000
A = [random.sample(range(m), m) for _ in range(m)]
B = [random.sample(range(m), m) for _ in range(m)]
match, proposals = gale_shapley(A, B)
rankA = [{x: i for i, x in enumerate(p)} for p in A]; rankB = [{x: i for i, x in enumerate(p)} for p in B]
husband = {r: p for p, r in match.items()}
blocking = sum(1 for p in range(m) for r in A[p][:rankA[p][match[p]]] if rankB[r][p] < rankB[r][husband[r]])
print(f"\nGale-Shapley, {m:,} a side: {proposals:,} proposals (at most {m * m:,}), blocking pairs: {blocking}")
avg_prop = sum(rankA[p][match[p]] for p in range(m)) / m + 1
avg_recv = sum(rankB[r][husband[r]] for r in range(m)) / m + 1
print(f"   average rank of partner: proposers {avg_prop:.1f}, receivers {avg_recv:.1f} (1 = first choice)")
match2, _ = gale_shapley(B, A)                  # same preferences, roles swapped
print(f"   roles swapped: the A side now averages {sum(rankA[a][b] for b, a in match2.items()) / m + 1:.1f}")

# 3. Hall's condition: a perfect matching exists iff every set S of people likes at least |S| jobs
likes = {"ann": {"db", "web"}, "bob": {"db"}, "cat": {"db", "web"}, "dan": {"ops", "ml", "db"}}
def max_matching(likes):
    owner = {}
    def augment(p, seen):
        for j in likes[p]:
            if j not in seen:
                seen.add(j)
                if j not in owner or augment(owner[j], seen): owner[j] = p; return True
        return False
    return sum(augment(p, set()) for p in likes), owner
size, owner = max_matching(likes)
worst = min((s for r in range(1, 5) for s in combinations(likes, r)),
            key=lambda s: len(set().union(*(likes[p] for p in s))) - len(s))
print(f"\nmaximum matching: {size} of {len(likes)} people; Hall violation: {', '.join(sorted(worst))} like only "
      f"{', '.join(sorted(set().union(*(likes[p] for p in worst))))}")

# 4. Euler's formula V - E + F = 2 on planar graphs, and the edge bounds that rule out K5 and K3,3
g = 10                                          # a g x g grid of points: count its parts directly
gV = g * g
gE = sum(1 for r in range(g) for c in range(g) for dr, dc in ((0, 1), (1, 0)) if r + dr < g and c + dc < g)
gF = (g - 1) ** 2 + 1                           # the unit squares, plus the outside
for name, V, E, F in (("cube", 8, 12, 6), ("octahedron", 6, 12, 8), ("icosahedron", 12, 30, 20), (f"{g}x{g} grid", gV, gE, gF)):
    print(f"{name:12} V - E + F = {V} - {E} + {F} = {V - E + F}")
print(f"K5:   E = {len(list(combinations(range(5), 2)))} > 3V - 6 = {3 * 5 - 6}, so K5 is not planar")
print(f"K3,3: E = {3 * 3} > 2V - 4 = {2 * 6 - 4} (no triangles), so K3,3 is not planar")

crown graph on 16 vertices (bipartite): greedy uses 8 colours in one order, 2 in another
random graph, max degree 29: greedy uses 9 colours (bound 30); largest-degree-first: 8

Gale-Shapley, 1,000 a side: 8,133 proposals (at most 1,000,000), blocking pairs: 0
   average rank of partner: proposers 8.1, receivers 122.9 (1 = first choice)
   roles swapped: the A side now averages 154.7

maximum matching: 3 of 4 people; Hall violation: ann, bob, cat like only db, web
cube         V - E + F = 8 - 12 + 6 = 2
octahedron   V - E + F = 6 - 12 + 8 = 2
icosahedron  V - E + F = 12 - 30 + 20 = 2
10x10 grid   V - E + F = 100 - 180 + 82 = 2
K5:   E = 10 > 3V - 6 = 9, so K5 is not planar
K3,3: E = 9 > 2V - 4 = 8 (no triangles), so K3,3 is not planar

Gale–Shapley needed 8,133 proposals, against the worst-case bound of a million. With random preferences the expected number is about n ln n, roughly 6,900 for n = 1,000, which is the coupon collector's problem of chapter 13 in disguise. The stability check went through every pair a proposer prefers to their partner and found no blocking pair.

The sum, product and complement rules

The sum rule: if a set splits into disjoint parts, its size is the sum of their sizes. The product rule: the number of sequences $(a_1,\dots ,a_k)$ with $a_i$ drawn from a set of size $n_i$ is $n_1{n}_2\cdots n_k$. It extends to choices that depend on earlier ones, provided the number of options at each step doesn't. That's the generalised product rule, which counts ordered choices without repetition below. There are ${62}^8\approx 2.2\times {10}^{14}$ passwords of 8 characters drawn from letters and digits, because each position independently takes one of 62 values.

The complement rule counts what you want by subtracting what you don't. Passwords with at least one digit are hard to count directly, since the digit could be anywhere and there could be several. Passwords with no digit are easy, ${52}^8$, so the answer is ${62}^8-{52}^8$. Whenever a problem says "at least one", try the complement first.

Permutations, combinations and the division rule

There are $n!$ ways to order n distinct objects, and $P(n,k)=n!/(n-k)!$ ways to choose an ordered sequence of k of them. To count unordered selections, use the division rule: if every object you want corresponds to exactly d of the things you counted, divide by d. Each set of k elements corresponds to k! ordered sequences, so the number of k-element subsets of an n-element set is

The division rule handles repetition too. The 11 letters of MISSISSIPPI can be ordered in 11! ways if you pretend the repeated letters are distinguishable, and every distinct word arises from $4!4!2!$ of those orderings, since the four I's, four S's and two P's can be shuffled among themselves. So there are $11!/(4!4!2!)=\mathrm{34,650}$ distinct words. People seated around a round table, where only relative positions matter, can be arranged in $n!/n=(n-1)!$ ways, because each seating corresponds to n rotations.

Bijections

The bijection rule (chapter 4) says that two sets in one-to-one correspondence have the same size. It turns a hard count into an easy one. Subsets of an n-element set correspond to bit strings of length n, one bit per element, so there are $2^n$ of them.

Stars and bars counts the ways to write 10 as $x_1+x_2+x_3$ with nonnegative integers. Represent a solution as 10 stars and 2 bars in a row, $x_1$ stars, a bar, $x_2$ stars, a bar, $x_3$ stars. Every arrangement of 10 stars and 2 bars is a solution and vice versa, so there are $\left(\genfrac{}{}{0ex}{}{12}{2}\right)=66$ solutions. In general, k nonnegative integers summing to n can be chosen in $\left(\genfrac{}{}{0ex}{}{n+k-1}{k-1}\right)$ ways. Similarly, monotone paths across a 6 × 6 grid correspond to strings of six R's and six U's, so there are $\left(\genfrac{}{}{0ex}{}{12}{6}\right)=924$.

Binomial coefficients and Pascal's triangle

The name comes from the binomial theorem:

Expanding the product of n factors (x + y) means choosing x or y from each factor, and the term $x^k{y}^{n-k}$ arises once for each way of choosing which k factors contribute x. Setting x = y = 1 gives ${\sum }_k\left(\genfrac{}{}{0ex}{}{n}{k}\right)=2^n$, which is another proof that a set has $2^n$ subsets.

Pascal's rule, $\left(\genfrac{}{}{0ex}{}{n}{k}\right)=\left(\genfrac{}{}{0ex}{}{n-1}{k-1}\right)+\left(\genfrac{}{}{0ex}{}{n-1}{k}\right)$, has a combinatorial proof, an argument that counts one set in two ways. Fix one element of the n-element set. A k-subset either contains it, and then the other k − 1 elements come from the remaining n − 1, or it doesn't, and all k come from the remaining n − 1. Combinatorial proofs are often shorter than algebra and explain why an identity holds.

Every count, checked

The example checks each rule against brute-force enumeration: ordered and unordered choices, the MISSISSIPPI words, stars and bars, the complement rule on a small alphabet, every 5-card poker hand classified one by one, and Pascal's rule and grid paths.

from itertools import combinations, permutations, product
from collections import Counter
from math import comb, perm, factorial

# 1. The basic rules, checked against brute force
print(f"ordered choices of 3 from 10: {sum(1 for _ in permutations(range(10), 3))} = P(10,3) = {perm(10, 3)}")
print(f"unordered choices of 3 from 10: {sum(1 for _ in combinations(range(10), 3))} = C(10,3) = {comb(10, 3)}")

# 2. The division rule: arrangements of MISSISSIPPI, generated without duplicates
def arrangements(counts):
    if not any(counts.values()): yield ""; return
    for ch in counts:
        if counts[ch]:
            counts[ch] -= 1
            for rest in arrangements(counts): yield ch + rest
            counts[ch] += 1
n_arr = sum(1 for _ in arrangements(Counter("MISSISSIPPI")))
print(f"arrangements of MISSISSIPPI: {n_arr:,} = 11!/(1! 4! 4! 2!) = {factorial(11) // (factorial(4) ** 2 * factorial(2)):,}")

# 3. Stars and bars: nonnegative solutions of x1 + x2 + x3 = 10
sols = sum(1 for x in range(11) for y in range(11 - x))          # z = 10 - x - y is forced
print(f"solutions of x1 + x2 + x3 = 10 in nonnegative integers: {sols} = C(12,2) = {comb(12, 2)}")

# 4. The complement rule: passwords with at least one digit, brute force on a small alphabet
alpha = "abc12"
with_digit = sum(1 for w in product(alpha, repeat=4) if any(ch.isdigit() for ch in w))
print(f"4-character strings over 'abc12' with a digit: {with_digit} = 5^4 - 3^4 = {5 ** 4 - 3 ** 4}")
print(f"8-character passwords over 62 symbols: {62 ** 8:,}; with at least one digit: {62 ** 8 - 52 ** 8:,}")

# 5. Every 5-card poker hand, classified
deck = range(52)                                                   # rank = card // 4, suit = card % 4
kinds = Counter()
for hand in combinations(deck, 5):
    shape = sorted(Counter(c >> 2 for c in hand).values(), reverse=True)
    if shape[0] == 4: kinds["four of a kind"] += 1
    elif shape[:2] == [3, 2]: kinds["full house"] += 1
    elif len({c & 3 for c in hand}) == 1: kinds["flush (any)"] += 1
    elif shape[0] == 3: kinds["three of a kind"] += 1
    elif shape[:2] == [2, 2]: kinds["two pair"] += 1
    elif shape[0] == 2: kinds["one pair"] += 1
print(f"\nall {comb(52, 5):,} five-card hands, counted one by one, against the formulas:")
formulas = {"four of a kind": 13 * 48, "full house": 13 * comb(4, 3) * 12 * comb(4, 2),
            "flush (any)": 4 * comb(13, 5), "three of a kind": 13 * comb(4, 3) * comb(12, 2) * 4 ** 2,
            "two pair": comb(13, 2) * comb(4, 2) ** 2 * 11 * 4, "one pair": 13 * comb(4, 2) * comb(12, 3) * 4 ** 3}
for k, f in formulas.items():
    print(f"   {k:16} {kinds[k]:>9,}   formula {f:>9,}")

# 6. Pascal's rule and the binomial theorem
row = [comb(12, k) for k in range(13)]
print(f"\nrow 12 of Pascal's triangle sums to {sum(row):,} = 2^12; C(12,6) = C(11,5) + C(11,6) = "
      f"{comb(11, 5)} + {comb(11, 6)} = {comb(12, 6)}")
paths = [[1] * 7 for _ in range(7)]
for i in range(1, 7):
    for j in range(1, 7): paths[i][j] = paths[i - 1][j] + paths[i][j - 1]
print(f"monotone paths across a 6 x 6 grid: {paths[6][6]} = C(12,6) = {comb(12, 6)}")

ordered choices of 3 from 10: 720 = P(10,3) = 720
unordered choices of 3 from 10: 120 = C(10,3) = 120
arrangements of MISSISSIPPI: 34,650 = 11!/(1! 4! 4! 2!) = 34,650
solutions of x1 + x2 + x3 = 10 in nonnegative integers: 66 = C(12,2) = 66
4-character strings over 'abc12' with a digit: 544 = 5^4 - 3^4 = 544
8-character passwords over 62 symbols: 218,340,105,584,896; with at least one digit: 164,880,377,053,440

all 2,598,960 five-card hands, counted one by one, against the formulas:
   four of a kind         624   formula       624
   full house           3,744   formula     3,744
   flush (any)          5,148   formula     5,148
   three of a kind     54,912   formula    54,912
   two pair           123,552   formula   123,552
   one pair         1,098,240   formula 1,098,240

row 12 of Pascal's triangle sums to 4,096 = 2^12; C(12,6) = C(11,5) + C(11,6) = 462 + 462 = 924
monotone paths across a 6 x 6 grid: 924 = C(12,6) = 924

The poker table is a good test of the method, because each formula is a product-rule argument that's easy to get subtly wrong. A full house is a rank for the triple (13 ways), three of its four suits ($\left(\genfrac{}{}{0ex}{}{4}{3}\right)$), a different rank for the pair (12) and two of its suits ($\left(\genfrac{}{}{0ex}{}{4}{2}\right)$): 3,744. Two pair needs care. The two pair ranks form an unordered set, $\left(\genfrac{}{}{0ex}{}{13}{2}\right)$, not 13 × 12, which would count every hand twice. Then come the suits for each pair and a fifth card from the 11 remaining ranks. Counting all 2,598,960 hands takes the example a few seconds and confirms all six formulas exactly. The "flush" line includes the 40 straight flushes, which a full hand ranking would separate out.

The pigeonhole principle

If n + 1 objects go into n boxes, some box holds at least two. More generally, if N objects go into k boxes, some box holds at least $\lceil N/k\rceil$. In the language of chapter 4: there's no one-to-one function from a larger finite set to a smaller one.

The principle proves several facts engineers rely on. Every hash function has collisions: a 64-bit hash maps infinitely many inputs to $2^{64}$ outputs, so some outputs are shared. That's unavoidable, and the design question is only whether collisions are hard to find. No lossless compressor can shrink every input: there are $2^n$ strings of n bits but only $2^n-1$ strings of all shorter lengths combined, so if a compressor shrank every n-bit string, two would compress to the same output and couldn't both be recovered. A compressor that shrinks some inputs must grow others, and the example shows zlib turning 10,000 random bytes into 10,011 while crushing 10,000 bytes of repetitive text into 99. And with 367 people, two must share a birthday. Chapter 12 shows the far more surprising fact that 23 people give better-than-even odds.

Order in every sequence

Pigeonhole arguments often hide in the choice of boxes. Erdős–Szekeres theorem: every sequence of $n^2+1$ distinct numbers contains an increasing or a decreasing subsequence of length n + 1. Proof. Label each position i with a pair $(a_i,d_i)$: the lengths of the longest increasing and longest decreasing subsequences ending at position i. For two positions i < j, if the j-th number is larger, then $a_j>a_i$, since it extends any increasing subsequence ending at i; if smaller, $d_j>d_i$. So no two positions have the same label. If every $a_i$ and $d_i$ were at most n, there would be only $n^2$ possible labels for $n^2+1$ positions, so two would coincide. Contradiction. ∎ The example confirms the theorem on 10,000 random sequences of 10 = 3² + 1 numbers, and shows a sequence of 9 that escapes it, so the bound is tight.

Inclusion–exclusion

The sum rule needs disjoint sets. For overlapping sets, adding the sizes counts the overlap twice, so subtract it: $|A\cup B|=|A|+|B|-|A\cap B|$. With three sets the pairwise overlaps are subtracted, but that removes the triple overlap once too often, so it's added back:

The general formula alternates over all intersections of 1, 2, 3, … of the sets. It's correct because an element lying in exactly m of the sets is counted $\left(\genfrac{}{}{0ex}{}{m}{1}\right)-\left(\genfrac{}{}{0ex}{}{m}{2}\right)+\left(\genfrac{}{}{0ex}{}{m}{3}\right)-\cdots =1$ times, by the binomial theorem applied to $(1-1{)}^m=0$. Inclusion–exclusion usually pairs with the complement rule. It's easier to count the objects that violate at least one condition and subtract them from the total than to count those that satisfy every condition directly.

Four applications

Divisibility. How many numbers from 1 to 1000 are divisible by none of 2, 3 and 5? Let $A_2,A_3,A_5$ be the multiples of each. Then $|A_2|=500$, $|A_3|=333$, $|A_5|=200$, the pairwise intersections are the multiples of 6, 10 and 15 (166, 100 and 66), and the triple intersection is the multiples of 30 (33). So the answer is 1000 − 1033 + 332 − 33 = 266. The same reasoning gives Euler's totient from chapter 6 as a product over the distinct primes dividing n:

Onto functions. Assign 6 distinct tasks to 3 workers so that nobody is idle. Of the $3^6$ assignments, subtract those that miss some particular worker, $3\cdot 2^6$, and add back those that miss two particular workers, $3\cdot 1^6$: 729 − 192 + 3 = 540.

Derangements. n people check their hats and get them back at random. How likely is it that nobody gets their own? A permutation with no fixed point is a derangement. Let $A_i$ be the permutations that fix person i. Any j particular people can be fixed in $(n-j)!$ ways, so inclusion–exclusion gives the number of derangements

The sum is the start of the series for $e^{-1}$, so $D_n/n!\to 1/e\approx 0.368$, and in fact $D_n$ is $n!/e$ rounded to the nearest integer. Remarkably, the probability barely depends on n: about 37% whether 9 people check hats or 9,000.

Password policies. How many 8-character passwords over letters and digits contain at least one lowercase letter, one uppercase letter and one digit? Count the violators. Passwords with no lowercase letter number ${36}^8$, with no uppercase ${36}^8$, with no digit ${52}^8$. Those lacking two kinds at once number ${10}^8$, ${26}^8$ and ${26}^8$, and none lack all three. So the policy admits ${62}^8-(2\cdot {36}^8+{52}^8)+({10}^8+2\cdot {26}^8)$ passwords, about 73% of all of them. A policy meant to strengthen passwords actually shrinks the space attackers must search, by about a quarter.

Checked by enumeration

The example verifies the compression count, Erdős–Szekeres, and each inclusion–exclusion formula above against brute force, including every permutation of up to 9 hats.

import math, random, zlib
from itertools import permutations, product

random.seed(10)
# 1. Pigeonhole: no lossless compressor can shrink every input
n = 16
print(f"bit strings of length {n}: {2 ** n:,}; strings of every shorter length combined: {2 ** n - 1:,}")
data = random.randbytes(10_000)
text = (b"the quick brown fox jumps over the lazy dog. " * 223)[:10_000]
print(f"zlib on 10,000 random bytes: {len(zlib.compress(data, 9)):,} bytes; on 10,000 bytes of English-like text: "
      f"{len(zlib.compress(text, 9)):,} bytes")

# 2. Erdos-Szekeres: any sequence of n^2 + 1 distinct numbers has a monotone run of length n + 1
def longest(seq, up):
    best = []
    for i, x in enumerate(seq):
        best.append(1 + max((best[j] for j in range(i) if (seq[j] < x) == up), default=0))
    return max(best)
trials = [random.sample(range(100), 10) for _ in range(10_000)]
print(f"10,000 random sequences of 10 = 3^2 + 1 numbers: every one has a monotone subsequence of length 4: "
      f"{all(max(longest(s, True), longest(s, False)) >= 4 for s in trials)}")
s9 = [7, 8, 9, 4, 5, 6, 1, 2, 3]
print(f"but 9 numbers can avoid it: {s9} has longest increasing {longest(s9, True)}, decreasing {longest(s9, False)}")

# 3. Inclusion-exclusion: numbers up to 1000 divisible by none of 2, 3, 5
brute = sum(1 for k in range(1, 1001) if k % 2 and k % 3 and k % 5)
ie = 1000 - (500 + 333 + 200) + (166 + 100 + 66) - 33
print(f"\n1..1000 divisible by none of 2, 3, 5: {brute} by brute force, {ie} by inclusion-exclusion")
phi = sum(1 for k in range(1, 361) if math.gcd(k, 360) == 1)
print(f"phi(360) = {phi} by counting, 360(1 - 1/2)(1 - 1/3)(1 - 1/5) = {360 * 1 * 2 * 4 // (2 * 3 * 5)}")

# 4. Onto functions: 6 tasks to 3 workers, nobody idle
onto = sum(1 for f in product(range(3), repeat=6) if len(set(f)) == 3)
print(f"assignments of 6 tasks to 3 workers using all 3: {onto} = 3^6 - 3 x 2^6 + 3 x 1^6 = {3**6 - 3 * 2**6 + 3}")

# 5. Derangements: nobody gets their own hat back
print("\n n   derangements   n!/e rounded   share of all n! orders")
for k in range(1, 10):
    d = sum(1 for p in permutations(range(k)) if all(p[i] != i for i in range(k)))
    print(f"{k:2}   {d:>12,}   {round(math.factorial(k) / math.e):>12,}   {d / math.factorial(k):.6f}")
print(f"1/e = {1 / math.e:.6f}")

# 6. Passwords of 8 characters with at least one lowercase, one uppercase and one digit
total = 62 ** 8
good = total - (36 ** 8 + 36 ** 8 + 52 ** 8) + (10 ** 8 + 26 ** 8 + 26 ** 8) - 0
small = sum(1 for w in product("aAb1", repeat=5) if any(c.islower() for c in w)
            and any(c.isupper() for c in w) and any(c.isdigit() for c in w))
print(f"\n8-character passwords with all three kinds: {good:,} of {total:,} ({good / total:.1%})")
print(f"same rule, alphabet 'aAb1', length 5: {small} by brute force, "
      f"{4**5 - (2**5 + 3**5 + 3**5) + (1**5 + 2**5 + 1**5)} by inclusion-exclusion")

bit strings of length 16: 65,536; strings of every shorter length combined: 65,535
zlib on 10,000 random bytes: 10,011 bytes; on 10,000 bytes of English-like text: 99 bytes
10,000 random sequences of 10 = 3^2 + 1 numbers: every one has a monotone subsequence of length 4: True
but 9 numbers can avoid it: [7, 8, 9, 4, 5, 6, 1, 2, 3] has longest increasing 3, decreasing 3

1..1000 divisible by none of 2, 3, 5: 266 by brute force, 266 by inclusion-exclusion
phi(360) = 96 by counting, 360(1 - 1/2)(1 - 1/3)(1 - 1/5) = 96
assignments of 6 tasks to 3 workers using all 3: 540 = 3^6 - 3 x 2^6 + 3 x 1^6 = 540

 n   derangements   n!/e rounded   share of all n! orders
 1              0              0   0.000000
 2              1              1   0.500000
 3              2              2   0.333333
 4              9              9   0.375000
 5             44             44   0.366667
 6            265            265   0.368056
 7          1,854          1,854   0.367857
 8         14,833         14,833   0.367882
 9        133,496        133,496   0.367879
1/e = 0.367879

8-character passwords with all three kinds: 159,655,911,367,680 of 218,340,105,584,896 (73.1%)
same rule, alphabet 'aAb1', length 5: 540 by brute force, 540 by inclusion-exclusion

The derangement table shows how fast the convergence is. The share of permutations with no fixed point is already 0.367882 at n = 8, against 1/e = 0.367879, because the error of a truncated alternating series is smaller than its first omitted term, here 1/9!. The same series explains why the rounded value $n!/e$ is always exactly right.

Sums that keep appearing

An arithmetic sum, such as the comparisons of insertion sort on reversed input, has the closed form ${\sum }_{i=1}^{n}i=n(n+1)/2$, which chapter 3 proved by induction. A geometric sum, with ratio r ≠ 1, is

Multiply the sum by r, subtract the original, and almost everything cancels. When r > 1 the largest term dominates. That's why an array that doubles its capacity copies fewer than twice as many elements in total as its final size: 1 + 2 + 4 + … + n < 2n. When r < 1 the first term dominates, and the infinite sum converges to 1/(1 − r).

The harmonic numbers $H_n=1+\frac{1}{2}+\cdots +\frac{1}{n}$ appear in the expected cost of randomised quicksort, the coupon collector's problem and the analysis of hashing. Comparing the sum with the integral of 1/x, the area under a staircase versus the area under a curve, gives $\ln (n+1)\le H_n\le 1+\ln n$, and more precisely

where $\gamma$ is Euler's constant. Factorials have Stirling's approximation, $n!\approx \sqrt{2\pi n}(n/e{)}^n$, with relative error about 1/(12n). Taking logarithms gives $\log n!=\mathrm{\Theta }(n\log n)$, the bound behind every comparison sort's lower limit in The Engineer's Workbench.

Asymptotic notation, precisely

f(n) = O(g(n)) means f grows no faster than g up to a constant: there are c and $n_0$ with $f(n)\le cg(n)$ for all $n\ge n_0$. $\mathrm{\Omega }$ is the matching lower bound, and $\mathrm{\Theta }$ means both. Little-o, f = o(g), means $f/g\to 0$. A practical test: if $\lim f(n)/g(n)$ exists and is a positive constant, then $f=\mathrm{\Theta }(g)$. Three habits avoid most mistakes. O is only an upper bound, so "this algorithm is O(n²)" is true of a linear algorithm too, and when you mean a tight bound, say $\mathrm{\Theta }$. The "=" in f = O(g) is not symmetric, since O(g) is really a set of functions. And the constants hidden by O decide which algorithm wins at the sizes you actually have.

From code to recurrence

A recursive algorithm's cost satisfies the same recursion as the code. Merge sort on n items sorts two halves and merges them with at most n comparisons, so its comparison count satisfies

To solve a recurrence, draw its recursion tree: one node per call, labelled with the work done outside the recursive calls.

The tree suggests T(n) = n log₂ n for powers of 2, and a guess can always be proved by induction, which is the substitution method (exercise 11.2). For n = 1,024 the example computes T exactly as 10,240 = n log₂ n. For other n, the floors and ceilings change the answer slightly but not the growth: 9,976 against 9,966 for n = 1,000.

The master theorem

Most divide-and-conquer algorithms split a problem of size n into a subproblems of size n/b and do $\mathrm{\Theta }(n^c)$ extra work. For

the master theorem compares c with ${\log }_{b}a$, the exponent that counts the leaves of the recursion tree:

The recursion tree explains all three cases. Level j has $a^j$ nodes each doing $(n/b^j{)}^c$ work, so the level totals form a geometric series with ratio $a/b^c$. If the ratio is below 1 the first term, the root, dominates. If it's exactly 1, all ${\log }_{b}n$ levels contribute equally. If it's above 1 the last level, the leaves, dominates, and there are $a^{{\log }_{b}n}=n^{{\log }_{b}a}$ of them.

Binary search (a = 1, b = 2, c = 0) is the middle case, with $\mathrm{\Theta }(\log n)$. Merge sort (2, 2, 1) is also the middle case, $\mathrm{\Theta }(n\log n)$. Karatsuba's multiplication replaces four half-size products with three (3, 2, 1), and the leaves dominate: $\mathrm{\Theta }(n^{{\log }_{2}3})=\mathrm{\Theta }(n^{1.585})$, beating the schoolbook $n^2$. Strassen's matrix multiplication uses seven half-size products instead of eight, giving $\mathrm{\Theta }(n^{{\log }_{2}7})=\mathrm{\Theta }(n^{2.807})$. The theorem doesn't apply when the extra work isn't a polynomial, such as $n/\log n$, or when the subproblems have unequal sizes. Then the Akra–Bazzi method or a direct look at the tree is needed.

Linear recurrences

Recurrences like $F_n=F_{n-1}+F_{n-2}$ have exact solutions. Guess $F_n=x^n$ and substitute, which gives the characteristic equation $x^2=x+1$, with roots $\phi =(1+\sqrt{5})/2\approx 1.618$ and $\psi =(1-\sqrt{5})/2\approx -0.618$. Every solution is a combination $A{\phi }^n+B{\psi }^n$, and the starting values $F_0=0,F_1=1$ fix the constants, giving Binet's formula:

Since $|\psi |<1$, the second term shrinks to nothing and $F_n$ is simply ${\phi }^n/\sqrt{5}$ rounded. That's mathematically exact for every n, but the example shows it failing in double-precision floating point at n = 71, where $F_{71}$ has 15 digits and the 53-bit significand can no longer carry the error below one half. For exact large Fibonacci numbers, use integer arithmetic with repeated squaring of a 2 × 2 matrix (chapter 17). The Towers of Hanoi recurrence $M_n=2M_{n-1}+1$ is solved by adding 1 to both sides: $M_n+1=2(M_{n-1}+1)$, so $M_n=2^n-1$.

Measured

The example checks the harmonic and Stirling approximations, computes the merge sort recurrence exactly, measures $T(n)$ divided by the master theorem's prediction for six recurrences, finds where Binet's formula fails in floating point, and counts Hanoi moves by recursion.

import math
from functools import lru_cache

# 1. Sums: harmonic numbers and Stirling's approximation
for n in (10, 1000, 10 ** 6):
    H = math.fsum(1 / k for k in range(1, n + 1))
    print(f"H({n:>9,}) = {H:9.5f};  H - ln n = {H - math.log(n):.5f}")
print(f"Euler's constant gamma = 0.57722 (the limit)")
for n in (10, 100, 1000):
    exact = math.lgamma(n + 1)                     # ln(n!)
    stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)
    print(f"n = {n:4}: n! / Stirling = {math.exp(exact - stirling):.6f}   (1 + 1/(12n) = {1 + 1 / (12 * n):.6f})")

# 2. Merge sort: T(n) = T(ceil(n/2)) + T(floor(n/2)) + n, T(1) = 0
@lru_cache(None)
def T(n): return 0 if n <= 1 else T((n + 1) // 2) + T(n // 2) + n
print("\nmerge sort recurrence: " + ", ".join(f"T({n:,}) = {T(n):,} vs n log2 n = {n * math.log2(n):,.0f}"
                                            for n in (1024, 1000, 10 ** 6)))

# 3. The master theorem: T(n) = a T(n/b) + n^c, measured at n = 2^k
def solve(a, b, c, n):
    return 1 if n <= 1 else a * solve(a, b, c, n // b) + n ** c
cases = [("binary search", 1, 2, 0, lambda n: math.log2(n) + 1),
         ("merge sort", 2, 2, 1, lambda n: n * math.log2(n)),
         ("Karatsuba", 3, 2, 1, lambda n: n ** math.log2(3)),
         ("4 halves, linear work", 4, 2, 1, lambda n: n ** 2),
         ("2 halves, quadratic work", 2, 2, 2, lambda n: n ** 2),
         ("Strassen", 7, 2, 2, lambda n: n ** math.log2(7))]
print(f"\n{'recurrence':26} {'predicted':>14}   T(n) / predicted at n = 2^10, 2^15, 2^20")
for name, a, b, c, f in cases:
    ratios = [solve(a, b, c, 2 ** k) / f(2 ** k) for k in (10, 15, 20)]
    e = math.log(a, b)                         # compare the leaves' exponent log_b a with c
    pred = ("log n" if (a, c) == (1, 0) else "n log n" if a == b ** c else
            (f"n^{e:.3f}" if abs(e - round(e)) > 1e-9 else f"n^{round(e)}") if a > b ** c else f"n^{c}")
    print(f"{name:26} {pred:>14}   " + "  ".join(f"{r:6.3f}" for r in ratios))

# 4. A linear recurrence: Fibonacci, and the closed form with floating point
phi = (1 + math.sqrt(5)) / 2
fib = [0, 1]
while len(fib) < 100: fib.append(fib[-1] + fib[-2])
bad = next(n for n in range(100) if round(phi ** n / math.sqrt(5)) != fib[n])
print(f"\nF(n) = round(phi^n / sqrt 5) holds exactly for n = 0..{bad - 1}; at n = {bad} floating point gives "
      f"{round(phi ** bad / math.sqrt(5)):,} instead of {fib[bad]:,}")

# 5. Towers of Hanoi: M(n) = 2 M(n-1) + 1
def hanoi(n, a="A", b="B", c="C"):
    return 0 if n == 0 else hanoi(n - 1, a, c, b) + 1 + hanoi(n - 1, c, b, a)
print("Hanoi moves: " + ", ".join(f"{n} discs {hanoi(n):,}" for n in (3, 10, 20)) + f"   (2^n - 1: {2**20 - 1:,} for 20)")

H(       10) =   2.92897;  H - ln n = 0.62638
H(    1,000) =   7.48547;  H - ln n = 0.57772
H(1,000,000) =  14.39273;  H - ln n = 0.57722
Euler's constant gamma = 0.57722 (the limit)
n =   10: n! / Stirling = 1.008365   (1 + 1/(12n) = 1.008333)
n =  100: n! / Stirling = 1.000834   (1 + 1/(12n) = 1.000833)
n = 1000: n! / Stirling = 1.000083   (1 + 1/(12n) = 1.000083)

merge sort recurrence: T(1,024) = 10,240 vs n log2 n = 10,240, T(1,000) = 9,976 vs n log2 n = 9,966, T(1,000,000) = 19,951,424 vs n log2 n = 19,931,569

recurrence                      predicted   T(n) / predicted at n = 2^10, 2^15, 2^20
binary search                       log n    1.000   1.000   1.000
merge sort                        n log n    1.100   1.067   1.050
Karatsuba                         n^1.585    2.965   2.995   2.999
4 halves, linear work                 n^2    1.999   2.000   2.000
2 halves, quadratic work              n^2    1.999   2.000   2.000
Strassen                          n^2.807    2.328   2.333   2.333

F(n) = round(phi^n / sqrt 5) holds exactly for n = 0..70; at n = 71 floating point gives 308,061,521,170,130 instead of 308,061,521,170,129
Hanoi moves: 3 discs 7, 10 discs 1,023, 20 discs 1,048,575   (2^n - 1: 1,048,575 for 20)

In the master theorem table, every ratio settles to a constant as n grows. That's exactly what $\mathrm{\Theta }$ promises, and the constants differ: about 3 for Karatsuba, 2.33 for Strassen, 1 in the limit for merge sort. The theorem says nothing about those constants, which is why Strassen's algorithm only beats ordinary matrix multiplication for matrices well into the hundreds or thousands of rows, depending on the implementation.

Sample spaces and the tree method

A sample space is the set of possible outcomes of an experiment, each with a probability, all summing to 1. An event is a set of outcomes, and its probability is the sum of theirs. When the experiment happens in stages, draw a tree: each branch is labelled with the probability of that step given everything before it, and an outcome's probability is the product along its path. The 6.042J text calls this the four-step method: find the outcomes, define the event, assign outcome probabilities, add them up. It sounds pedantic, and it settles arguments that intuition can't.

Take the Monty Hall problem. A car is behind one of three doors and goats are behind the others. You pick door 1. The host, who knows where the car is, opens a different door that has a goat, and offers you the chance to switch to the remaining door. Should you?

The tree settles it: switching wins with probability 2/3. The common intuition, "two doors left, so 50–50", ignores the fact that the host's choice depends on where the car is. The example confirms 2/3 exactly by walking the tree with fractions, and approximately by playing 100,000 games.

Conditional probability

The probability of A given that B has happened is

Knowing B shrinks the sample space to B, and the formula rescales so that B's outcomes sum to 1. Rearranged, it's the product rule, $\mathrm{Pr}[A\cap B]=\mathrm{Pr}[B]\mathrm{Pr}[A\mid B]$, which is what multiplying along a tree branch does. If events $B_1,\dots ,B_k$ split the sample space, the law of total probability adds up the branches: $\mathrm{Pr}[A]={\sum }_i\mathrm{Pr}[B_i]\mathrm{Pr}[A\mid B_i]$.

What you condition on must be exactly what you learned, including how you learned it. With two children, each equally likely to be a boy or a girl, $\mathrm{Pr}[\text{two boys}\mid \text{at least one boy}]=1/3$, since three equally likely families remain: BB, BG, GB. But $\mathrm{Pr}[\text{two boys}\mid \text{the older is a boy}]=1/2$. The information sounds similar, and it leaves different sample spaces.

Bayes' rule

Writing $\mathrm{Pr}[A\cap B]$ in two ways gives Bayes' rule:

It converts what you can measure, how often the evidence appears when the cause is present, into what you want, how likely the cause is given the evidence. A diagnostic test detects a condition 95% of the time and gives a false positive 5% of the time. One person in a hundred has the condition. If you test positive, the chance that you have it is

Only about 16%. The test is accurate, but the condition is rare, so the healthy 99% produce far more false positives than the sick 1% produce true ones. Ignoring that is the base rate fallacy. A second independent positive test raises the probability to about 79%, because the first result's posterior becomes the second one's prior. Monitoring has the same structure. An alert that is 99% accurate, watching for a failure that occurs in 0.1% of minutes, will page you mostly for nothing, and the fix is better specificity or a combination of signals, not a lower threshold.

Independence

A and B are independent if $\mathrm{Pr}[A\cap B]=\mathrm{Pr}[A]\mathrm{Pr}[B]$, equivalently if learning B doesn't change the probability of A. Independence is a numerical fact to check, not a matter of whether events "seem related". Roll two dice. "The first die is even" is independent of "the sum is 7", since whatever the first die shows, exactly one value of the second makes 7. It is not independent of "the sum is 8". Given an even first die, the second die can always complete 8 (6, 4 or 2), so the sum is 8 with probability 1/6, against 5/36 overall. The example checks both by counting. Assuming independence when it doesn't hold is one of the costliest mistakes in reliability engineering. Disks from the same batch fail together, and replicas in the same rack lose power together, so "two failures at once" is far likelier than the product of the individual probabilities.

Birthdays

How many people must be in a room before a shared birthday is more likely than not? Counting the complement, with 365 equally likely days, n people all have different birthdays with probability

and this drops below 1/2 at n = 23. The surprise comes from counting pairs, not people. 23 people form 253 pairs, and each pair has a 1-in-365 chance of matching. Chapter 14 applies the same approximation to random identifiers and hash values, where it decides how many bits you need.

Computed and simulated

The example computes each result exactly with fractions, then simulates it: Monty Hall, the test with its base rate across a million simulated people, the birthday problem, the two-child puzzle, and independence for two pairs of dice events.

import random
from fractions import Fraction as F
from itertools import product

random.seed(12)
# 1. Monty Hall, exactly: walk every branch of the tree with its probability
win_stay = win_switch = F(0)
for car in (1, 2, 3):                          # you always pick door 1
    goats = [d for d in (2, 3) if d != car]    # the host opens a goat door you didn't pick
    for opened in goats:
        p = F(1, 3) * F(1, len(goats))
        other = next(d for d in (1, 2, 3) if d not in (1, opened))
        win_stay += p * (car == 1); win_switch += p * (car == other)
print(f"Monty Hall exactly: stay wins {win_stay}, switch wins {win_switch}")
N = 100_000; stay = switch = 0
for _ in range(N):
    car, pick = random.randint(1, 3), random.randint(1, 3)
    opened = random.choice([d for d in (1, 2, 3) if d not in (car, pick)])
    other = next(d for d in (1, 2, 3) if d not in (pick, opened))
    stay += car == pick; switch += car == other
print(f"Monty Hall simulated ({N:,} games): stay {stay / N:.3f}, switch {switch / N:.3f}")

# 2. Bayes' rule and the base rate: a 95%-accurate test for a condition 1 in 100 people have
prior, sens, fpr = F(1, 100), F(95, 100), F(5, 100)
post = sens * prior / (sens * prior + fpr * (1 - prior))
post2 = sens * post / (sens * post + fpr * (1 - post))
print(f"\nP(condition | positive test) = {float(post):.3f}; after a second independent positive: {float(post2):.3f}")
people = [(random.random() < 0.01) for _ in range(1_000_000)]
pos = [sick for sick in people if random.random() < (0.95 if sick else 0.05)]
print(f"simulated million people: {len(pos):,} positives, of whom {sum(pos) / len(pos):.3f} have the condition")

# 3. The birthday problem
def p_shared(n):
    p_distinct = F(1)
    for i in range(n): p_distinct *= F(365 - i, 365)
    return 1 - p_distinct
print("\nP(some shared birthday): " + ", ".join(f"{n} people {float(p_shared(n)):.3f}" for n in (10, 23, 30, 50, 70)))
hits = sum(len(set(random.randrange(365) for _ in range(23))) < 23 for _ in range(N))
print(f"simulated rooms of 23 people: {hits / N:.3f} contain a shared birthday")

# 4. Conditioning on what you're told: two children
families = list(product("BG", repeat=2))       # (older, younger), equally likely
at_least_one = [f for f in families if "B" in f]
older_boy = [f for f in families if f[0] == "B"]
print(f"\nP(two boys | at least one boy) = {F(sum(f == ('B', 'B') for f in at_least_one), len(at_least_one))}; "
      f"P(two boys | older is a boy) = {F(sum(f == ('B', 'B') for f in older_boy), len(older_boy))}")

# 5. Independence is a numerical fact, not an intuition: two dice
dice = list(product(range(1, 7), repeat=2))
def P(event): return F(sum(1 for d in dice if event(d)), len(dice))
A = lambda d: d[0] % 2 == 0
for s in (7, 8):
    B = lambda d, s=s: sum(d) == s
    both = P(lambda d: A(d) and B(d))
    print(f"first die even, sum = {s}: P(both) = {both}, P(A)P(B) = {P(A) * P(B)} -> "
          f"{'independent' if both == P(A) * P(B) else 'dependent'}")

Monty Hall exactly: stay wins 1/3, switch wins 2/3
Monty Hall simulated (100,000 games): stay 0.332, switch 0.668

P(condition | positive test) = 0.161; after a second independent positive: 0.785
simulated million people: 58,709 positives, of whom 0.160 have the condition

P(some shared birthday): 10 people 0.117, 23 people 0.507, 30 people 0.706, 50 people 0.970, 70 people 0.999
simulated rooms of 23 people: 0.510 contain a shared birthday

P(two boys | at least one boy) = 1/3; P(two boys | older is a boy) = 1/2
first die even, sum = 7: P(both) = 1/12, P(A)P(B) = 1/12 -> independent
first die even, sum = 8: P(both) = 1/12, P(A)P(B) = 5/72 -> dependent

The simulations land close to the exact values, within the sampling error chapter 19 quantifies: 0.668 against 2/3, 0.160 against 0.161, 0.510 against 0.507. That's the right role for simulation. It checks a calculation, or explores a problem too messy to calculate. It doesn't replace the calculation when one is available, because it can't tell you whether 0.668 means 2/3 or 0.67.

Random variables and expectation

A random variable X assigns a number to each outcome of an experiment. Its expectation is the probability-weighted average of its values:

A fair die has expectation (1 + 2 + … + 6)/6 = 7/2. The expectation needn't be a possible value: no roll shows 3.5. An indicator variable $I_A$ is 1 if event A happens and 0 otherwise, so $𝔼[I_A]=\mathrm{Pr}[A]$. Indicators are how you turn questions about how many things happen into questions about how likely each one is.

Linearity of expectation

For any random variables X and Y, and constants a and b,

Proof. Write both sides as sums over outcomes $\omega$: ${\sum }_{\omega }(aX(\omega )+bY(\omega ))\mathrm{Pr}[\omega ]$ splits into $a{\sum }_{\omega }X(\omega )\mathrm{Pr}[\omega ]+b{\sum }_{\omega }Y(\omega )\mathrm{Pr}[\omega ]$. ∎ The proof never mentions independence, and it isn't needed. That's what makes linearity so powerful: write a complicated count as a sum of indicators, which may depend on each other in tangled ways, and add up their probabilities.

Fixed points. n people get their hats back in a random order. How many get their own, on average? Let $I_i$ indicate that person i does. Each person gets their own hat with probability 1/n, so the expected number is $n\cdot 1/n=1$, for every n. The indicators are dependent, since if n − 1 people have their own hats the last one must too, and it doesn't matter.

Inversions. How many swaps does insertion sort make on a random permutation? One per inversion, a pair of positions in the wrong order. Each of the $\left(\genfrac{}{}{0ex}{}{n}{2}\right)$ pairs is inverted with probability 1/2, so the expected number is $n(n-1)/4$. That makes insertion sort quadratic on average, not just in the worst case.

Running maxima. Scanning an array to find its maximum, how often does the "best so far" change? Position i holds a new maximum exactly when its value is the largest of the first i, which for a random order has probability 1/i. The expected number of updates is $1+\frac{1}{2}+\cdots +\frac{1}{n}=H_n\approx \ln n$, only about 7.5 for a thousand values. That's also the expected number of times a "keep the best candidate" loop has to write.

Empty bins. Throw n balls into n bins at random. Each bin stays empty with probability $(1-1/n{)}^n\approx 1/e$, so about 37% of the bins are empty, and the same fraction of hash buckets is unused when a table has as many keys as buckets.

Waiting times and the coupon collector

If each trial succeeds independently with probability p, the number of trials up to and including the first success has a geometric distribution with expectation 1/p. One trial always happens, and with probability 1 − p you're back where you started, so $𝔼=1+(1-p)𝔼$, which solves to 1/p.

The coupon collector's problem asks how many random draws it takes to collect all n kinds of coupon. Split the process into phases. Once you have i kinds, each draw is new with probability (n − i)/n, so that phase takes n/(n − i) draws on average. By linearity the total expectation is

Variance

The expectation says nothing about spread. The variance measures it as the expected squared distance from the mean, and the standard deviation $\sigma$ is its square root:

A die has $𝔼[X^2]=91/6$, so $\mathrm{Var}[X]=91/6-49/4=35/12$. Unlike expectation, variance adds only for independent variables, or more generally uncorrelated ones: $\mathrm{Var}[X+Y]=\mathrm{Var}[X]+\mathrm{Var}[Y]$. So the sum of 10 dice has variance 350/12 ≈ 29.2 and standard deviation about 5.4, not 10 times the single die's 1.7. The spread of a sum grows like the square root of the number of terms, which is why averages are stable (chapter 19). For dependent variables, add twice their covariance.

Expectations hide the spread, and in engineering the spread is often what matters. The coupon collector's average for 50 kinds is 225 draws, but in the example one run needed 706. Request latency has the same character: a service with a fine mean latency can have a tail that ruins the experience for one request in a hundred, which is why production systems measure percentiles, not just averages (The Engineer's Watch).

Simulated

The example estimates each expectation above by repeating the experiment thousands of times and compares the average with the formula.

import random
from fractions import Fraction as F
from statistics import mean, pvariance

random.seed(13)
def H(n): return sum(1 / k for k in range(1, n + 1))

# 1. Fixed points of a random permutation: expected number 1, whatever n (linearity of expectation)
for n in (10, 100):
    counts = []
    for _ in range(20_000):
        p = list(range(n)); random.shuffle(p)
        counts.append(sum(p[i] == i for i in range(n)))
    print(f"fixed points of a random permutation of {n:3}: mean {mean(counts):.3f}, variance {pvariance(counts):.3f}")

# 2. Coupon collector: draws needed to see all 50 kinds
runs = []
for _ in range(10_000):
    seen, draws = set(), 0
    while len(seen) < 50: seen.add(random.randrange(50)); draws += 1
    runs.append(draws)
print(f"coupon collector, 50 kinds: mean {mean(runs):.1f} draws (n H_n = {50 * H(50):.1f}), "
      f"fewest {min(runs)}, most {max(runs)}")

# 3. How often does a running maximum change? Expected H_n times
updates = []
for _ in range(10_000):
    best, k = -1, 0
    for x in random.sample(range(10 ** 6), 1000):
        if x > best: best, k = x, k + 1
    updates.append(k)
print(f"updates to a running maximum over 1,000 random values: mean {mean(updates):.2f} (H_1000 = {H(1000):.2f})")

# 4. Inversions: the swaps insertion sort makes on a random permutation
n = 100
inv = []
for _ in range(2_000):
    p = random.sample(range(n), n)
    inv.append(sum(p[i] > p[j] for i in range(n) for j in range(i + 1, n)))
print(f"inversions in a random permutation of {n}: mean {mean(inv):.0f} (n(n-1)/4 = {n * (n - 1) // 4})")

# 5. Empty bins: n balls into n bins leave about n/e bins empty
n = 1000
empties = [n - len(set(random.randrange(n) for _ in range(n))) for _ in range(2_000)]
print(f"{n} balls into {n} bins: mean empty bins {mean(empties):.1f} (n(1 - 1/n)^n = {n * (1 - 1 / n) ** n:.1f})")

# 6. Variance: exact for one die, additive for independent dice
die = range(1, 7)
EX = F(sum(die), 6); VX = F(sum(x * x for x in die), 6) - EX ** 2
sums = [sum(random.randint(1, 6) for _ in range(10)) for _ in range(50_000)]
print(f"one die: E = {EX}, Var = {VX} = {float(VX):.3f};  sum of 10 dice: mean {mean(sums):.2f} (35), "
      f"variance {pvariance(sums):.2f} ({float(10 * VX):.2f})")

fixed points of a random permutation of  10: mean 0.999, variance 1.012
fixed points of a random permutation of 100: mean 0.995, variance 0.999
coupon collector, 50 kinds: mean 224.1 draws (n H_n = 225.0), fewest 96, most 706
updates to a running maximum over 1,000 random values: mean 7.48 (H_1000 = 7.49)
inversions in a random permutation of 100: mean 2478 (n(n-1)/4 = 2475)
1000 balls into 1000 bins: mean empty bins 368.0 (n(1 - 1/n)^n = 367.7)
one die: E = 7/2, Var = 35/12 = 2.917;  sum of 10 dice: mean 35.02 (35), variance 29.23 (29.17)

Every simulated mean lands within a fraction of a percent of its formula, and the fixed-point experiment adds something the formula didn't promise: the variance is about 1 as well. That's a hint of the Poisson distribution of chapter 14, which the number of fixed points approaches as n grows.

A field guide to distributions

The Poisson distribution is the binomial's limit when n is large, p is small and $\lambda =np$ stays fixed: $\mathrm{Pr}[X=k]=e^{-\lambda }{\lambda }^k/k!$. It models the number of keys in one bucket of a hash table, or requests arriving in one second from many independent users. The geometric distribution is memoryless: having failed 100 times doesn't change the chance of success on the next try. The normal distribution appears whenever many independent contributions add up, a fact chapter 19 makes precise with the central limit theorem. Be wary of using it in the far tails, where it can be wrong by orders of magnitude.

Markov and Chebyshev

Markov's inequality. For a random variable X ≥ 0 and any a > 0,

Proof. $𝔼[X]\ge a\cdot \mathrm{Pr}[X\ge a]$, since X contributes at least a whenever X ≥ a and never contributes a negative amount. ∎ It's weak, but it needs nothing except the mean and non-negativity. If an algorithm's expected running time is 1 second, at most 10% of runs take 10 seconds or more. Restarting runs that exceed twice the expected time therefore gives a guarantee, since each attempt succeeds with probability at least 1/2.

Chebyshev's inequality. For any X with mean $\mu$ and variance ${\sigma }^2$,

Apply Markov to the non-negative variable $(X-\mu {)}^2$ with $a=t^2$. It says a variable is rarely more than a few standard deviations from its mean: at most 1/9 of the probability lies beyond 3σ, whatever the distribution. Applied to an average of n independent samples, whose variance is ${\sigma }^2/n$, it proves the law of large numbers: averages converge to the mean.

Chernoff and Hoeffding

When X is a sum of independent pieces, its tails are far thinner than Chebyshev admits. They're exponentially small. Hoeffding's inequality: if X is a sum of n independent variables each between 0 and 1, then