The Engineer's Workbench

Why contiguous memory wins

An array stores its elements next to each other, so the address of element i is the start address plus i times the element size: constant-time access with no pointers to follow. Contiguity also suits the hardware from stage 04. Reading one element pulls its neighbours into the same cache line, so a loop over an array streams through memory at full speed, while a loop over a linked list jumps to a new, probably uncached, address at every step. Python's list is an array, of references rather than values, and so are Java's ArrayList, C++'s vector and Go's slices.

Growing without paying for it

A fixed-size block can't grow in place, because the memory after it may belong to something else. So when a dynamic array is full, it allocates a new, larger block, copies every element across, and frees the old one. A single append can therefore cost as much as the whole array. What saves it is how much larger the new block is.

class DynamicArray:
    """A growable array on top of fixed-size storage: when it's full, allocate more and copy."""
    def __init__(self, grow=lambda cap: cap * 2, shrink_at=4):
        self.cap, self.n, self.slots = 1, 0, [None]
        self.grow, self.shrink_at, self.copies, self.resizes = grow, shrink_at, 0, 0
    def _resize(self, new_cap):
        slots = [None] * new_cap
        for i in range(self.n):                    # the expensive part: copy every element
            slots[i] = self.slots[i]
            self.copies += 1
        self.slots, self.cap, self.resizes = slots, new_cap, self.resizes + 1
    def append(self, x):
        if self.n == self.cap:
            self._resize(max(self.cap + 1, self.grow(self.cap)))
        self.slots[self.n] = x
        self.n += 1
    def pop(self):
        self.n -= 1
        x, self.slots[self.n] = self.slots[self.n], None
        if self.cap > 1 and self.n <= self.cap // self.shrink_at:
            self._resize(max(1, self.cap // 2))    # give memory back, but not too eagerly
        return x
    def __getitem__(self, i):
        if not 0 <= i < self.n:
            raise IndexError(i)
        return self.slots[i]

N = 50_000
print(f"appending {N:,} items:")
for label, grow in (("double (x2)", lambda c: c * 2), ("grow by half (x1.5)", lambda c: c + c // 2),
                    ("add 1,000 slots", lambda c: c + 1000)):
    a = DynamicArray(grow)
    for i in range(N):
        a.append(i)
    print(f"   {label:20} {a.resizes:4} resizes, {a.copies:>10,} copies = {a.copies / N:6.2f} per append")

print("pushing and popping at the boundary, 1,000 times:")
for label, shrink_at in (("shrink when half full", 2), ("shrink when a quarter full", 4)):
    a = DynamicArray(shrink_at=shrink_at)
    for i in range(1024):
        a.append(i)
    a.copies = a.resizes = 0
    for _ in range(1000):
        a.append(0); a.pop()
    print(f"   {label:27} {a.resizes:5} resizes, {a.copies:>9,} copies")

class RingBuffer:
    """A fixed-capacity queue: head and tail indices chase each other around one array."""
    def __init__(self, capacity):
        self.slots, self.head, self.size = [None] * capacity, 0, 0
    def push(self, x):
        if self.size == len(self.slots):
            raise OverflowError("full")
        self.slots[(self.head + self.size) % len(self.slots)] = x
        self.size += 1
    def pop(self):
        x = self.slots[self.head]
        self.head, self.size = (self.head + 1) % len(self.slots), self.size - 1
        return x

q, moved, out_ring, out_list, lst = RingBuffer(1000), 0, [], [], []
for i in range(1000):
    q.push(i); lst.append(i)
for i in range(1000, 6000):                        # a queue under steady load: one in, one out
    out_ring.append(q.pop()); q.push(i)
    moved += len(lst) - 1                          # list.pop(0) shifts every remaining element left
    out_list.append(lst.pop(0)); lst.append(i)
print(f"queue of 1,000 under 5,000 push/pop pairs: ring buffer moved 0 elements, "
      f"list.pop(0) moved {moved:,}; same order out: {out_ring == out_list}")

appending 50,000 items:
   double (x2)            16 resizes,     65,535 copies =   1.31 per append
   grow by half (x1.5)    27 resizes,    122,904 copies =   2.46 per append
   add 1,000 slots        50 resizes,  1,225,050 copies =  24.50 per append
pushing and popping at the boundary, 1,000 times:
   shrink when half full        2000 resizes, 2,048,000 copies
   shrink when a quarter full      1 resizes,     1,024 copies
queue of 1,000 under 5,000 push/pop pairs: ring buffer moved 0 elements, list.pop(0) moved 4,995,000; same order out: True

With doubling, 50,000 appends caused 16 resizes and 65,535 copies, 1.31 copies per append. The copies form a geometric series, 1 + 2 + 4 + ... + 32,768, which always totals less than twice the final size. That's the amortised cost: some appends are expensive, but averaged over any sequence of appends, each one costs a constant. Growing by half is just as constant, at 2.46 copies per append, and leaves less memory unused. Growing by a fixed 1,000 slots is a different kind of thing. It resizes every 1,000 appends and copies everything each time, so the total is quadratic: 1.2 million copies here, 24.5 per append, and the per-append figure keeps rising as the array grows.

An accounting argument gives the same result without the series. Charge every append three units: one pays for writing the element, and two go into a savings account. When the array doubles from k to 2k, the k/2 elements appended since the last resize have saved k units, exactly enough to pay for copying all k elements. The account never goes negative, so three units per append always covers the work.

Real libraries pick different factors for different trade-offs. Java's ArrayList grows by half; GNU's C++ library doubles a vector, while Microsoft's grows it by half; Go doubles small slices and grows large ones more gently; CPython over-allocates a list by about an eighth plus a small constant. Smaller factors waste less memory, larger ones copy less, and every one of them is a multiple, never a fixed step.

Shrinking without thrashing

An array that never shrinks wastes memory after a burst of deletions, so dynamic arrays usually give space back. The obvious rule, halve the capacity when the array falls to half full, has a trap the second block of output shows. At the boundary, one append fills the array and forces a doubling, and the next pop empties it to half and forces a halving. Alternating the two copies the whole array every time: 2,000 resizes and over two million copies for 1,000 pairs of operations. Shrinking only when the array falls to a quarter full, and then only to half, leaves room to move in both directions. The same 2,000 operations then cost one resize. This gap between the grow threshold and the shrink threshold, called hysteresis, appears wherever a system resizes itself, from hash tables to autoscaling groups.

Queues: the ring buffer

A queue adds at one end and removes at the other, and a plain array is bad at removing from the front: list.pop(0) shifts every remaining element left by one. Under a steady load of 1,000 queued items, 5,000 removals moved nearly five million elements. A ring buffer fixes this by never moving anything. It keeps a head index and a size, writes new items at (head + size) mod capacity, and advances the head on removal, wrapping around to the start of the array when it reaches the end. Both operations are constant time, with no copying, and the memory is allocated once. That predictability is why ring buffers sit inside network cards, audio pipelines, logging systems and the Linux kernel's io_uring interface. When a queue must also grow, use a ring buffer that doubles (exercise 1.3), or Python's collections.deque, a linked list of fixed-size blocks that is constant time at both ends.

From key to slot

A hash function turns a key into a large number, and the number modulo the table size picks a slot. If every key landed in its own slot, a lookup would be one array access. They don't. Collisions arrive much sooner than intuition suggests, by the same arithmetic as the birthday paradox: with 365 slots, 23 keys are more likely than not to include a collision (exercise 2.1). So the real design question is what to do when two keys want the same slot, and there are two families of answers.

Separate chaining makes each slot a list, and a lookup scans the list for its key. It degrades gracefully: at an average of one key per slot, chains are short. Open addressing keeps every key in the array itself. On a collision, it probes other slots in a fixed sequence, the next one for linear probing, until it finds the key or an empty slot. Probing avoids the pointers and allocations of chains and reads consecutive memory, which modern hardware rewards, but it is far more sensitive to how full the table is.

import random

MASK64 = (1 << 64) - 1
def mix(k):
    """A deterministic integer hash (splitmix64's finaliser): nearby keys land far apart."""
    k = (k + 0x9E3779B97F4A7C15) & MASK64
    k = ((k ^ (k >> 30)) * 0xBF58476D1CE4E5B9) & MASK64
    k = ((k ^ (k >> 27)) * 0x94D049BB133111EB) & MASK64
    return k ^ (k >> 31)

class ChainedMap:
    def __init__(self, hash_fn=mix, buckets=8):
        self.h, self.buckets, self.n, self.compares = hash_fn, [[] for _ in range(buckets)], 0, 0
    def _bucket(self, key):
        return self.buckets[self.h(key) % len(self.buckets)]
    def put(self, key, value):
        b = self._bucket(key)
        for pair in b:
            if pair[0] == key:
                pair[1] = value; return
        b.append([key, value]); self.n += 1
        if self.n > 0.75 * len(self.buckets):                  # keep chains short: double and rehash
            old, self.buckets = self.buckets, [[] for _ in range(2 * len(self.buckets))]
            for chain in old:
                for k, v in chain:
                    self._bucket(k).append([k, v])
    def get(self, key):
        for k, v in self._bucket(key):
            self.compares += 1
            if k == key: return v
        return None

class ProbingMap:
    """Open addressing with linear probing: on a collision, try the next slot."""
    EMPTY, GONE = object(), object()
    def __init__(self, capacity):
        self.slots, self.probes = [self.EMPTY] * capacity, 0
    def put(self, key, value):
        i = mix(key) % len(self.slots)
        while self.slots[i] is not self.EMPTY and self.slots[i] is not self.GONE and self.slots[i][0] != key:
            i = (i + 1) % len(self.slots)
        self.slots[i] = (key, value)
    def get(self, key):
        i = mix(key) % len(self.slots)
        while True:
            self.probes += 1
            s = self.slots[i]
            if s is self.EMPTY: return None
            if s is not self.GONE and s[0] == key: return s[1]
            i = (i + 1) % len(self.slots)
    def delete(self, key):
        i = mix(key) % len(self.slots)
        while self.slots[i] is not self.EMPTY:
            if self.slots[i] is not self.GONE and self.slots[i][0] == key:
                self.slots[i] = self.GONE; return            # a tombstone keeps later probe chains intact
            i = (i + 1) % len(self.slots)

random.seed(2)
CAP = 20_000
print("linear probing: average probes per lookup (Knuth's formulas in brackets)")
for alpha in (0.25, 0.5, 0.75, 0.9):
    m, keys = ProbingMap(CAP), random.sample(range(10**9), int(CAP * alpha))
    for k in keys: m.put(k, k)
    m.probes = 0
    for k in keys: m.get(k)
    hit = m.probes / len(keys)
    m.probes = 0
    for k in random.sample(range(10**9, 2 * 10**9), 5000): m.get(k)
    miss = m.probes / 5000
    print(f"   load {alpha:4}: hit {hit:5.2f} [{(1 + 1 / (1 - alpha)) / 2:5.2f}]   "
          f"miss {miss:6.2f} [{(1 + 1 / (1 - alpha) ** 2) / 2:6.2f}]")

m, keys = ProbingMap(CAP), random.sample(range(10**9), CAP // 2)
for k in keys: m.put(k, k)
for k in keys[: CAP // 4]: m.delete(k)                        # delete half; the slots become tombstones
m.probes = 0
for k in random.sample(range(10**9, 2 * 10**9), 5000): m.get(k)
print(f"after deleting half the keys, a miss still costs {m.probes / 5000:.2f} probes: "
      "tombstones count as occupied until a rebuild")

chains = {}
for name, fn in (("good hash (mix)", mix), ("weak hash (k mod 1024)", lambda k: k % 1024)):
    hm = ChainedMap(fn, buckets=1024)
    hm.buckets = [[] for _ in range(1024)]
    attack = [i * 1024 for i in range(2000)]                   # keys chosen so k mod 1024 == 0
    for k in attack:
        b = hm.buckets[fn(k) % 1024]; b.append([k, k])         # insert without resizing, to show the shape
    for k in attack: hm.get(k)
    chains[name] = max(len(b) for b in hm.buckets)
    print(f"{name:24} 2,000 crafted keys: longest chain {chains[name]:5,}, "
          f"{hm.compares / len(attack):7.1f} comparisons per lookup")

linear probing: average probes per lookup (Knuth's formulas in brackets)
   load 0.25: hit  1.16 [ 1.17]   miss   1.40 [  1.39]
   load  0.5: hit  1.51 [ 1.50]   miss   2.50 [  2.50]
   load 0.75: hit  2.54 [ 2.50]   miss   8.63 [  8.50]
   load  0.9: hit  5.35 [ 5.50]   miss  47.56 [ 50.50]
after deleting half the keys, a miss still costs 2.53 probes: tombstones count as occupied until a rebuild
good hash (mix)          2,000 crafted keys: longest chain     8,     2.0 comparisons per lookup
weak hash (k mod 1024)   2,000 crafted keys: longest chain 2,000,  1000.5 comparisons per lookup

Load factor is everything

The first block measures linear probing at four load factors, the fraction of slots in use, against the formulas Donald Knuth derived in 1963. They agree closely. A successful lookup stays cheap for a long time: 1.16 probes at a quarter full, 2.54 at three-quarters. Unsuccessful lookups, which must probe until they hit an empty slot, explode: 2.5 probes at half full, 8.6 at three-quarters, and 47.6 at 90%. The cause is clustering. Keys that collide fill neighbouring slots, so the run of full slots grows, and every key that hashes anywhere into the run extends it further. That's why tables resize long before they're full: Java's HashMap at 0.75, CPython's dict at two-thirds. Resizing allocates a table twice the size and reinserts every key, because each key's slot depends on the table size. As with chapter 1's array, doubling makes the cost amortised constant.

Deletion leaves tombstones

Deleting from an open-addressed table can't simply empty the slot. A key placed further along the probe sequence, because this slot was occupied when it arrived, would become unreachable, since a lookup stops at the first empty slot it meets. So deletion leaves a tombstone: a marker that lookups step over and inserts may reuse. The example's second line shows the cost. After deleting half the keys, the table holds as many live keys as a quarter-full one, but a miss still costs 2.53 probes, the half-full figure, because tombstones count as occupied. They are cleared only when the table is rebuilt. Alternative schemes avoid tombstones by shifting later entries back on deletion, or by Robin Hood hashing, which keeps probe lengths even.

Hash flooding

Every guarantee above assumes keys are spread evenly. The last two lines show what happens when an adversary chooses them. With a weak hash, the key modulo the table size, 2,000 keys that are all multiples of 1,024 land in one bucket, and every lookup scans a chain of 2,000: 1,000.5 comparisons where the mixing hash needs 2.0. This isn't theoretical. In December 2011, researchers showed that a single crafted HTTP request with thousands of colliding form field names could keep a web server's CPU busy for minutes, across PHP, Java, Python, Ruby and others, because those languages parsed request parameters into hash tables with predictable hashes. The fixes are why modern hash tables look the way they do. Python 3.4 adopted SipHash, a keyed hash with a random key per process (PEP 456), so attackers can't predict collisions. Java 8's HashMap turns any long chain into a balanced tree (chapter 4), capping the damage at logarithmic.

What the library versions add

Production hash tables combine these ideas with engineering for the hardware. CPython's dict keeps entries in a compact array in insertion order and a separate small index of slots, which is why dictionaries remember insertion order (guaranteed since Python 3.7). Google's Swiss tables and Rust's hashbrown store a byte of each hash in a metadata array and check sixteen slots at once with SIMD instructions. All of them still depend on the same three conditions: a good hash, a bounded load factor, and keys an attacker can't aim.

A tree hiding in an array

A binary heap is a binary tree with one rule, the heap property: every parent is no larger than its children (for a min-heap). It doesn't keep the whole tree sorted, only each parent above its children, which is much cheaper to maintain and still puts the smallest item at the root. The second trick is that the tree is always complete, filled level by level from the left, so it fits in an array with no gaps and no pointers. The root is at index 0, and the children of index i are at 2i + 1 and 2i + 2.

Both operations restore the heap property along one path, so both cost the height of the tree, about log₂ n. To push, append the item at the end and sift it up, swapping it with its parent while it's smaller. To pop the minimum, take the root, move the last item into its place, and sift it down, swapping it with its smaller child while either child is smaller.

import random

class MinHeap:
    """A binary heap stored in a list: the children of index i live at 2i+1 and 2i+2."""
    def __init__(self, items=()):
        self.a, self.compares = list(items), 0
        for i in range(len(self.a) // 2 - 1, -1, -1):      # bottom-up heapify: O(n), not O(n log n)
            self._down(i)
    def _less(self, i, j):
        self.compares += 1
        return self.a[i] < self.a[j]
    def _up(self, i):
        while i > 0 and self._less(i, (i - 1) // 2):
            p = (i - 1) // 2
            self.a[i], self.a[p] = self.a[p], self.a[i]
            i = p
    def _down(self, i):
        n = len(self.a)
        while True:
            l, r, small = 2 * i + 1, 2 * i + 2, i
            if l < n and self._less(l, small): small = l
            if r < n and self._less(r, small): small = r
            if small == i: return
            self.a[i], self.a[small] = self.a[small], self.a[i]
            i = small
    def push(self, x):
        self.a.append(x); self._up(len(self.a) - 1)
    def pop(self):
        last = self.a.pop()
        if not self.a: return last
        top, self.a[0] = self.a[0], last
        self._down(0)
        return top
    def __len__(self):
        return len(self.a)

random.seed(3)
N = 100_000
print(f"building a heap of {N:,} items, counting comparisons:")
for label, data in (("random order", [random.random() for _ in range(N)]),
                    ("descending order", list(range(N, 0, -1)))):       # every push climbs to the root
    one_by_one = MinHeap()
    for x in data: one_by_one.push(x)
    bottom_up = MinHeap(data)
    print(f"   {label:17} one push at a time {one_by_one.compares:>9,}   bottom-up heapify {bottom_up.compares:>9,}")
print("heapsort matches sorted():", [bottom_up.pop() for _ in range(len(bottom_up))] == sorted(data))

# The same heap drives shortest-path search. A* adds a guess of the remaining distance.
W = H = 60
def search(heuristic, walls):
    start, goal = (0, 0), (W - 1, H - 1)
    best, frontier, expanded = {start: 0}, MinHeap([(heuristic(start), 0, start)]), 0
    while frontier:
        _, neg_g, (x, y) = frontier.pop()
        g = -neg_g
        if g > best.get((x, y), float("inf")): continue     # a stale entry: a shorter way was found
        expanded += 1
        if (x, y) == goal: return g, expanded
        for nx, ny in ((x + 1, y), (x - 1, y), (x, y + 1), (x, y - 1)):
            if 0 <= nx < W and 0 <= ny < H and (nx, ny) not in walls and g + 1 < best.get((nx, ny), float("inf")):
                best[(nx, ny)] = g + 1
                frontier.push((g + 1 + heuristic((nx, ny)), -(g + 1), (nx, ny)))   # ties go to the deeper cell
    return None, expanded

def manhattan(p): return (W - 1 - p[0]) + (H - 1 - p[1])  # never overestimates on a 4-connected grid
for seed in range(100):                                    # the first maze, in a fixed order, with a way through
    random.seed(seed)
    maze = {(x, y) for x in range(W) for y in range(H) if random.random() < 0.3} - {(0, 0), (W - 1, H - 1)}
    if search(lambda p: 0, maze)[0] is not None:
        break
print(f"shortest path across a {W}x{H} grid, cells expanded:")
for label, walls in (("open ground", set()), ("maze, 30% walls", maze)):
    (d_len, d_exp), (a_len, a_exp) = search(lambda p: 0, walls), search(manhattan, walls)
    print(f"   {label:16} Dijkstra {d_exp:5,} cells, A* {a_exp:5,} cells, same length: {d_len == a_len} ({a_len} steps)")

building a heap of 100,000 items, counting comparisons:
   random order      one push at a time   228,554   bottom-up heapify   188,339
   descending order  one push at a time 1,468,946   bottom-up heapify   199,978
heapsort matches sorted(): True
shortest path across a 60x60 grid, cells expanded:
   open ground      Dijkstra 3,600 cells, A*   119 cells, same length: True (118 steps)
   maze, 30% walls  Dijkstra 2,465 cells, A* 2,005 cells, same length: True (128 steps)

Building a heap in linear time

The first block compares two ways to turn 100,000 items into a heap. Pushing them one at a time costs up to log n per item, and on descending input, where every new item is the smallest yet and climbs all the way to the root, it costs 1.47 million comparisons. On random input it's far cheaper, because most new items stop near the bottom. Bottom-up heapify works the other way: it sifts down every internal node, from the last one back to the root. It costs about 200,000 comparisons on either input, because it's guaranteed linear. Half the nodes are leaves and never move, a quarter can sift down at most one level, an eighth at most two, and the sum of those distances is less than n. When you have all the items up front, heapify them; Python's heapq.heapify does exactly this.

Popping everything from a heap returns the items in order, which is heapsort: O(n log n) in the worst case, in place, with no extra memory. It's rarely the library sort anyway. It isn't stable, and its jumps between parents and children defeat the cache, so in practice it serves as a safety net: introsort, the usual C++ std::sort, switches to heapsort when quicksort starts going badly (chapter 5).

Priority queues in the wild

Anything that repeatedly asks for the most urgent item is a priority queue. Operating system schedulers pick the next task, event simulations process the earliest event, timers fire the soonest deadline, and "the top 100 of ten million" keeps a heap of 100 (exercise 3.2). One operation the basic heap lacks is changing the priority of an item already inside it, which Dijkstra's algorithm wants. Production code either keeps an index from items to positions (an indexed heap), or does what the example does: push the item again with its new priority and skip the stale copy when it's popped. Python's heapq has no decrease-key at all, so that second approach is standard there.

A*: Dijkstra with a sense of direction

Dijkstra's algorithm, which Sorted Again's graphs module teaches, expands nodes in order of their distance from the start, g. It finds the shortest path, but it searches outwards in every direction equally, including away from the goal. A* orders the queue by g + h instead, where h is an estimate of the distance still to go. If the estimate never overestimates, it's called admissible, and A* is guaranteed to find a shortest path while expanding far fewer nodes. Dijkstra is simply A* with h = 0.

The second block shows how much depends on the estimate. On open ground, the Manhattan distance is exactly right, and A* expands only the 119 cells of the path itself, where Dijkstra expands all 3,600. In a maze with 30% walls, the estimate is optimistic, since walls force detours it can't see, so A* saves less: 2,005 cells against 2,465. Both always find the same length. Two details make it work: the heap holds (g + h, -g, cell), so among equally promising cells the deeper one goes first, and a popped entry that's worse than a distance already found is skipped as stale. Route planners for real road networks go further: stronger estimates built from precomputed distances to landmarks, and preprocessing schemes such as contraction hierarchies that answer continental queries thousands of times faster than plain Dijkstra.

What a search tree is for

Sorted Again's trees module covers binary search trees: every key in a node's left subtree is smaller than it, every key in the right subtree is larger, an in-order traversal visits the keys in sorted order, and a search follows a single path from the root. Hash tables are faster for exact lookups, so the reason to use a tree is everything a hash table can't do: iterate in order, find every key in a range, find the next key above or below a value, or the k-th smallest. The catch is that every operation costs the height of the tree, and nothing in a plain BST controls the height.

import random

class BST:
    """A plain binary search tree: fine on random input, a linked list on sorted input."""
    def __init__(self):
        self.root = None                                   # nodes are [key, left, right]
    def insert(self, key):
        if self.root is None:
            self.root = [key, None, None]; return
        node = self.root
        while True:                                        # iterative: a degenerate tree is very deep
            side = 1 if key < node[0] else 2
            if node[side] is None:
                node[side] = [key, None, None]; return
            node = node[side]
    def height(self):
        h, level = 0, [self.root] if self.root else []
        while level:
            h, level = h + 1, [c for n in level for c in (n[1], n[2]) if c]
        return h

class AVL:
    """Keeps every node's two subtree heights within one of each other, by rotating."""
    class Node:
        __slots__ = ("key", "left", "right", "h")
        def __init__(self, key): self.key, self.left, self.right, self.h = key, None, None, 1
    def __init__(self):
        self.root, self.rotations = None, 0
    @staticmethod
    def _h(n): return n.h if n else 0
    def _fix(self, n): n.h = 1 + max(self._h(n.left), self._h(n.right))
    def _rotate_right(self, y):
        x = y.left
        y.left, x.right = x.right, y
        self._fix(y); self._fix(x); self.rotations += 1
        return x
    def _rotate_left(self, x):
        y = x.right
        x.right, y.left = y.left, x
        self._fix(x); self._fix(y); self.rotations += 1
        return y
    def _insert(self, n, key):
        if n is None: return self.Node(key)
        if key < n.key: n.left = self._insert(n.left, key)
        else: n.right = self._insert(n.right, key)
        self._fix(n)
        balance = self._h(n.left) - self._h(n.right)
        if balance > 1:                                    # left-heavy
            if key >= n.left.key: n.left = self._rotate_left(n.left)    # left-right case
            return self._rotate_right(n)
        if balance < -1:                                   # right-heavy
            if key < n.right.key: n.right = self._rotate_right(n.right) # right-left case
            return self._rotate_left(n)
        return n
    def insert(self, key):
        self.root = self._insert(self.root, key)
    def inorder(self):
        out, stack, n = [], [], self.root
        while stack or n:
            while n: stack.append(n); n = n.left
            n = stack.pop(); out.append(n.key); n = n.right
        return out

N = 4095                                                   # 2^12 - 1: a perfect tree would have height 12
random.seed(4)
heights = {}
for label, keys in (("sorted keys", list(range(N))), ("random keys", random.sample(range(10**6), N))):
    plain, avl = BST(), AVL()
    for k in keys:
        plain.insert(k); avl.insert(k)
    heights[label] = (plain.height(), avl.root.h)
    print(f"{label:12} {N:,} inserts: plain BST height {plain.height():5}   "
          f"AVL height {avl.root.h:3} ({avl.rotations:,} rotations)   AVL in order: {avl.inorder() == sorted(keys)}")
worst_plain, worst_avl = heights["sorted keys"]
print(f"a search walks one root-to-leaf path, so on sorted input it costs up to {worst_plain:,} steps "
      f"in the plain tree and {worst_avl} in the AVL tree")

sorted keys  4,095 inserts: plain BST height  4095   AVL height  12 (4,083 rotations)   AVL in order: True
random keys  4,095 inserts: plain BST height    24   AVL height  15 (2,767 rotations)   AVL in order: True
a search walks one root-to-leaf path, so on sorted input it costs up to 4,095 steps in the plain tree and 12 in the AVL tree

How a tree degenerates

Inserting 4,095 keys in sorted order into a plain BST gives a tree of height 4,095. Each new key is larger than everything before it, so it becomes the right child of the previous one, and the "tree" is a linked list with extra pointers. A search for the largest key walks all 4,095 nodes. Sorted input isn't an unlucky accident: auto-incrementing ids, timestamps, alphabetical imports and already-sorted files are some of the most common inputs there are. Random keys do much better, a height of 24 here, but a structure whose speed depends on the order the data happens to arrive in can't be trusted.

Rotations

A rotation is a local restructuring that changes which of two nodes is the parent while keeping the in-order sequence, and so the search-tree property, intact. It touches three pointers and runs in constant time.

AVL trees: never more than one apart

The AVL tree, the first self-balancing tree (Adelson-Velsky and Landis, 1962), keeps one invariant: at every node, the heights of the two subtrees differ by at most one. After an insertion, the code walks back up the path, updates heights, and wherever a node's balance reaches two, fixes it with a rotation. There are four cases. In the left-left and right-right cases, a single rotation fixes it. In the left-right and right-left cases, the inner child is rotated first to turn it into one of the simple cases. The invariant bounds the height at about 1.44 log₂ n, so every search, insert and delete is logarithmic in the worst case.

The output shows the payoff. On the same 4,095 sorted keys that flattened the plain tree, the AVL tree has height 12, exactly what a perfect tree of 4,095 nodes would have, after 4,083 rotations: sorted input triggers one at almost every insertion, and the rotations pack the keys perfectly. On random keys it's 15, well within its bound, and an in-order traversal returns the keys sorted in both cases, confirming the rotations never broke the search-tree property.

The other balanced trees

Red-black trees keep a looser invariant, colouring nodes red or black and allowing a height of up to 2 log₂(n + 1). In return, insertions and deletions need fewer rotations, which makes them the usual choice for general-purpose ordered maps: Java's TreeMap, most implementations of C++'s std::map, and many structures inside the Linux kernel. B-trees balance in a different way, with hundreds of keys per node to suit disk pages, which is why databases use them (The Engineer's Ledger, chapter 6). Skip lists get logarithmic expected time from random coin flips instead of rotations; Redis uses them for sorted sets, and LevelDB and RocksDB use them for their memtables. Python's standard library has no balanced tree at all. Where you need one, the popular sortedcontainers package uses a list of sorted lists, which performs very well in practice.

The floor under every comparison sort

A comparison sort learns about the order only by asking "is this item smaller than that one?" Picture every possible run of the algorithm as a tree of those questions, with a yes branch and a no branch at each one. Each leaf is a final ordering, and since any of the n! orderings of n distinct items might be the right one, the tree needs at least n! leaves. A binary tree with n! leaves has height at least log₂(n!), so some input must cost that many comparisons. By Stirling's approximation, log₂(n!) is about n log₂ n − 1.44n, which for 2,000 items is 19,053.

Seven sorts, four inputs

The example implements the classic algorithms and counts the comparisons each makes on 2,000 items in four arrangements: random, already sorted, reversed, and drawn from only ten distinct values.

import functools, math, random

C = 0
def less(a, b):
    global C
    C += 1
    return a < b

def insertion_sort(a):
    for i in range(1, len(a)):
        j = i
        while j > 0 and less(a[j], a[j - 1]):            # slide each item left into place
            a[j], a[j - 1] = a[j - 1], a[j]
            j -= 1

def merge_sort(a):
    if len(a) <= 1: return a
    mid = len(a) // 2
    left, right = merge_sort(a[:mid]), merge_sort(a[mid:])
    out, i, j = [], 0, 0
    while i < len(left) and j < len(right):
        if less(right[j], left[i]): out.append(right[j]); j += 1
        else: out.append(left[i]); i += 1                # ties take the left item: stable
    return out + left[i:] + right[j:]

def quick_sort(a, pivot_rule):
    stack = [(0, len(a) - 1)]                            # an explicit stack: no recursion limit
    while stack:
        lo, hi = stack.pop()
        if lo >= hi: continue
        p = pivot_rule(lo, hi)
        a[p], a[hi] = a[hi], a[p]
        store = lo
        for i in range(lo, hi):                          # partition around the pivot at a[hi]
            if less(a[i], a[hi]):
                a[i], a[store] = a[store], a[i]; store += 1
        a[store], a[hi] = a[hi], a[store]
        stack += [(lo, store - 1), (store + 1, hi)]

def quick_sort_3way(a, rng):
    stack = [(0, len(a) - 1)]
    while stack:
        lo, hi = stack.pop()
        if lo >= hi: continue
        pivot = a[rng.randint(lo, hi)]
        lt, i, gt = lo, lo, hi                           # < pivot | == pivot | unknown | > pivot
        while i <= gt:
            if less(a[i], pivot): a[lt], a[i] = a[i], a[lt]; lt += 1; i += 1
            elif less(pivot, a[i]): a[i], a[gt] = a[gt], a[i]; gt -= 1
            else: i += 1
        stack += [(lo, lt - 1), (gt + 1, hi)]            # the equal block is already in place

def heap_sort(a):
    n = len(a)
    def down(i, end):
        while True:
            l, r, big = 2 * i + 1, 2 * i + 2, i
            if l < end and less(a[big], a[l]): big = l
            if r < end and less(a[big], a[r]): big = r
            if big == i: return
            a[i], a[big] = a[big], a[i]; i = big
    for i in range(n // 2 - 1, -1, -1): down(i, n)       # max-heap
    for end in range(n - 1, 0, -1):
        a[0], a[end] = a[end], a[0]; down(0, end)        # move the max to the back

def timsort(a):
    def compare(x, y):                                   # one call is one comparison
        global C
        C += 1
        return (x > y) - (x < y)
    a.sort(key=functools.cmp_to_key(compare))

N = 2000
random.seed(5)
inputs = {"random": random.sample(range(10**6), N), "sorted": list(range(N)),
          "reversed": list(range(N, 0, -1)), "10 values": [random.randrange(10) for _ in range(N)]}
rng = random.Random(9)
algorithms = {
    "insertion": insertion_sort,
    "merge": lambda a: a.__setitem__(slice(None), merge_sort(a)),
    "quick, first pivot": lambda a: quick_sort(a, lambda lo, hi: lo),
    "quick, random pivot": lambda a: quick_sort(a, lambda lo, hi: rng.randint(lo, hi)),
    "quick, 3-way": lambda a: quick_sort_3way(a, rng),
    "heap": heap_sort,
    "Timsort (sorted())": timsort,
}
print(f"comparisons to sort {N:,} items  (lower bound for any comparison sort: log2({N}!) = {math.lgamma(N + 1) / math.log(2):,.0f})")
print(f"{'':22}" + "".join(f"{k:>12}" for k in inputs))
for name, sort in algorithms.items():
    row = []
    for data in inputs.values():
        a = list(data); C = 0
        sort(a)
        assert a == sorted(data)
        row.append(C)
    print(f"{name:22}" + "".join(f"{c:>12,}" for c in row))

class Rec:
    """A record compared by its key only, so records with equal keys are ties to the sort."""
    def __init__(self, key, pos): self.key, self.pos = key, pos
    def __lt__(self, other): return self.key < other.key
def keeps_order(recs):
    return all(a.pos < b.pos for a, b in zip(recs, recs[1:]) if a.key == b.key)

recs = [Rec(random.randrange(5), i) for i in range(30)]
by_merge = merge_sort(list(recs))
by_heap = list(recs); heap_sort(by_heap)
print(f"\nequal keys keep their original order?  merge: {keeps_order(by_merge)}   "
      f"heap: {keeps_order(by_heap)}   sorted(): {keeps_order(sorted(recs))}")

comparisons to sort 2,000 items  (lower bound for any comparison sort: log2(2000!) = 19,053)
                            random      sorted    reversed   10 values
insertion                1,014,248       1,999   1,999,000     883,754
merge                       19,422      10,864      11,088      18,803
quick, first pivot          23,845   1,999,000   1,999,000     205,652
quick, random pivot         23,182      27,573      24,797     208,800
quick, 3-way                37,460      39,061      38,206      13,057
heap                        37,789      39,159      35,964      35,302
Timsort (sorted())          19,303       1,999       1,999      13,822

equal keys keep their original order?  merge: True   heap: False   sorted(): True

Insertion sort slides each item left into place. On random input it's quadratic, a million comparisons, and on reversed input two million, but on sorted input it makes one comparison per item and stops: 1,999. That makes it the fastest sort for tiny or nearly sorted arrays, which is why every hybrid sort uses it for small pieces.

Merge sort splits the list in half, sorts each half, and merges them. It's close to the floor on random input, 19,422 against 19,053, never worse than n log n, and stable, but it needs a second array as large as the input. Its merging of sorted runs is also the basis of external sorting, where data is too large for memory.

Quicksort partitions around a pivot and recurses on each side, in place. Its weakness is the pivot. Always taking the first item is quadratic on sorted and reversed input, 1,999,000 comparisons, exactly n(n−1)/2, because every partition peels off one item. A random pivot makes that vanishingly unlikely, at around 23,000 to 27,000 comparisons on any order. But the same row shows a second weakness: with only ten distinct values, a two-way partition puts all the equal items on one side, and it climbs back to 208,800. Three-way partitioning groups the items equal to the pivot in the middle and never touches them again, which brings the duplicates case down to 13,057. It pays for that with two comparisons per item on distinct keys, about 37,500.

Heapsort (chapter 3) is the steadiest: about twice the floor on every input, in place, with an O(n log n) worst case. It's slow in practice because its jumps through the array defeat the cache, and it isn't stable.

Timsort, which Tim Peters designed for Python in 2002, is a merge sort that first looks for runs of data already in order, ascending or descending, then merges them carefully. On random input it lands near the floor, 19,303. On sorted and reversed input it finds one long run and finishes in n−1 comparisons, and it does well with duplicates too. It's stable, and it's the sort behind Python's sorted(), Java's sort for objects since Java 7, and JavaScript's Array.prototype.sort in V8.

Stability

A sort is stable if items that compare equal keep their original relative order. The last line tests it with records sorted by a key that has only five values: merge sort and Python's sorted() preserve the original order within each key, and heapsort scrambles it. Stability matters whenever you sort records by one field and care about the order from a previous sort: a table sorted by date, then by customer, shows each customer's orders still in date order only if the second sort is stable. It also makes multi-key sorting simple. Sort by the least important key first, then by the next, each with a stable sort (exercise 5.3), which is exactly how the radix sort in chapter 6 works.

What libraries actually do

No single algorithm wins every column, so standard libraries combine them. Introsort, the usual C++ std::sort, runs quicksort, switches to heapsort if the recursion gets suspiciously deep, and finishes small pieces with insertion sort: quicksort's speed, heapsort's worst case, insertion sort's small-array efficiency. Pattern-defeating quicksort, used by Go's sort package since version 1.19 and long used for Rust's unstable sort, adds detection of sorted and repetitive input. Stable library sorts are typically Timsort or a relative of it. The practical rule follows. Use the library sort unless you have measured a reason not to, know whether it's stable, and reach for a specialised algorithm only for a specific need: a partial sort for the top k, quickselect for a median, an external merge sort for data larger than memory.

Counting sort

If every key is an integer between 0 and k−1, you don't need to compare anything. Count how many items have each key, turn the counts into starting positions with a running total, and place each item at its key's next free position. That's two passes over the items and one over the counts: n + k steps, linear when k isn't much larger than n. Placing items in their original order makes it stable, a property the next section depends on entirely.

import random

def counting_sort(items, key, k):
    """Stable sort of items whose keys are integers in [0, k): n + k steps, no comparisons."""
    counts = [0] * k
    for x in items: counts[key(x)] += 1
    total = 0
    for v in range(k):                                    # prefix sums: where each key's run starts
        counts[v], total = total, total + counts[v]
    out = [None] * len(items)
    for x in items:                                       # walk forwards: equal keys keep their order
        out[counts[key(x)]] = x
        counts[key(x)] += 1
    return out, len(items) * 2 + k

def radix_sort(nums, digit_pass=counting_sort):
    """LSD radix sort of 32-bit integers, one byte at a time, least significant first."""
    work = 0
    for shift in (0, 8, 16, 24):
        nums, w = digit_pass(nums, lambda x: (x >> shift) & 0xFF, 256)
        work += w
    return nums, work

def unstable_pass(items, key, k):
    """The same buckets, filled back to front: still groups by digit, but reverses equal keys."""
    counts = [0] * k
    for x in items: counts[key(x)] += 1
    total = 0
    for v in range(k): counts[v], total = total, total + counts[v]
    ends = [counts[v + 1] if v + 1 < k else len(items) for v in range(k)]
    out = [None] * len(items)
    for x in items:
        ends[key(x)] -= 1; out[ends[key(x)]] = x
    return out, len(items) * 2 + k

C = 0
def merge_sort(a):
    global C
    if len(a) <= 1: return a
    l, r = merge_sort(a[:len(a) // 2]), merge_sort(a[len(a) // 2:])
    out, i, j = [], 0, 0
    while i < len(l) and j < len(r):
        C += 1
        if r[j] < l[i]: out.append(r[j]); j += 1
        else: out.append(l[i]); i += 1
    return out + l[i:] + r[j:]

random.seed(6)
N = 100_000
nums = [random.getrandbits(32) for _ in range(N)]
result, work = radix_sort(nums)
merged = merge_sort(nums)
print(f"sorting {N:,} random 32-bit integers:")
print(f"   merge sort       {C:>10,} comparisons")
print(f"   LSD radix sort   {work:>10,} steps (4 passes of n + 256), no comparisons; correct: {result == sorted(nums)}")

ages = [random.randrange(120) for _ in range(N)]
out, work = counting_sort(ages, key=lambda x: x, k=120)
print(f"   counting sort of {N:,} ages 0-119: {work:,} steps; correct: {out == sorted(ages)}")
print(f"   counting sort of 64-bit ids would need a count array of 2^64 entries: "
      f"{(2**64 * 8) / 2**60:,.0f} EiB of memory")

floats = [random.random() for _ in range(N)]
buckets = [[] for _ in range(N)]
for f in floats: buckets[int(f * N)].append(f)            # uniform data: about one item per bucket
longest = max(len(b) for b in buckets)
out = [f for b in buckets for f in sorted(b)]
print(f"   bucket sort of {N:,} uniform floats: longest bucket {longest}; correct: {out == sorted(floats)}")

broken, _ = radix_sort(nums, digit_pass=unstable_pass)
first_bad = next(i for i in range(N - 1) if broken[i] > broken[i + 1])
print(f"\nradix sort with an unstable digit pass: sorted? {broken == sorted(nums)} "
      f"(first out-of-order pair at position {first_bad:,})")

sorting 100,000 random 32-bit integers:
   merge sort        1,536,586 comparisons
   LSD radix sort      801,024 steps (4 passes of n + 256), no comparisons; correct: True
   counting sort of 100,000 ages 0-119: 200,120 steps; correct: True
   counting sort of 64-bit ids would need a count array of 2^64 entries: 128 EiB of memory
   bucket sort of 100,000 uniform floats: longest bucket 9; correct: True

radix sort with an unstable digit pass: sorted? False (first out-of-order pair at position 1)

Sorting 100,000 ages between 0 and 119 takes 200,120 steps, with no comparisons at all. The next line shows the limit. The count array has one slot per possible key, so sorting 64-bit ids directly would need 2⁶⁴ counters, 128 exbibytes of memory. Counting sort is for small key ranges: ages, grades, status codes, bytes.

Radix sort

Radix sort extends counting sort to large keys by sorting them one digit at a time, using a stable counting sort for each digit. The least significant digit (LSD) version starts from the ones digit and works up. It looks backwards, but it works, and the diagram shows why.

Take 329 and 355. They share a hundreds digit, so the hundreds pass can't tell them apart: a stable pass leaves them in the order they arrived in. They arrived in the right order because the tens pass sorted them by their tens digits, 2 before 5. So the final order of two numbers that tie on a digit is decided by the earlier passes, and it survives only if every later pass is stable. That's the whole correctness argument for LSD radix sort.

The example sorts 100,000 random 32-bit integers in four passes of one byte each, 801,024 steps against merge sort's 1,536,586 comparisons, and the gap widens with n, since radix sort is linear and merge sort isn't. The last line is the warning. Swap in a digit pass that groups by digit correctly but reverses the order of equal digits, and the output is wrong from position 1 onwards, with no error raised. The bug is silent: every pass looks correct on its own.

Bucket sort

Bucket sort assumes the keys are spread evenly over a known range. It splits the range into n buckets, drops each item into its bucket in one pass, sorts each tiny bucket, and concatenates them. With 100,000 uniformly random floats and 100,000 buckets, the fullest bucket held 9 items, so the whole sort is linear in expectation. The assumption carries all the weight. If the data clusters, say most values near zero, most items land in a few buckets, and the sort degrades to whatever sorts those buckets.

When the floor doesn't apply

None of this contradicts chapter 5. The n log n bound counts comparisons, and these sorts don't compare: they use the structure of the keys, which a comparison sort isn't allowed to see. That's also why libraries rarely use them by default. A general-purpose sort can't know that your keys are small integers or uniformly spread. Where those facts are known, they win clearly: GPU libraries sort with radix sort, databases sort fixed-width keys with it, suffix-array construction depends on it, and counting sort is the inner loop of many of them.

The price of exact

Answering "have I seen this user before?" exactly means storing every user. Answering "how many distinct visitors today?" exactly means remembering every visitor until midnight. At a few hundred million items that's gigabytes per question, per machine. Hashing offers a way out. A good hash function makes any item look like a random number, and random numbers have statistics you can rely on. Each structure in this chapter stores only those statistics, in an amount of memory set by the error you'll accept rather than by how much data arrives.

Bloom filters: membership with one-sided error

A Bloom filter, invented by Burton Bloom in 1970, is an array of m bits, all zero, and k hash functions. To add an item, compute its k positions and set those bits to 1. To query an item, check the same k positions. If any of them is 0, the item was definitely never added, because adding it would have set that bit. If all of them are 1, the item was probably added, or other items happened to set all k bits between them, which is a false positive. There are never false negatives.

After n insertions, the chance that a given bit is still 0 is about e^−kn/m, so the false-positive rate is about (1 − e^−kn/m)^k. With too few hash functions each query checks too little; with too many the array fills with ones. The best choice is k = (m/n) ln 2, which leaves half the bits set. At that k, the bits you need depend only on the error rate, not on the size of the items: about 1.44 log₂(1/p) bits per key, which is 9.6 bits for 1% and 14.4 bits for 0.1%. Each tenfold improvement costs another 4.8 bits per key. Implementations don't compute k independent hashes. Adam Kirsch and Michael Mitzenmacher showed that two hashes combined as h₁ + i·h₂ work as well, which is what the example does with one blake2b digest split in half.

The classic use is guarding an expensive lookup that usually fails. LSM-tree databases such as RocksDB and Cassandra keep a Bloom filter for each file on disk, so a read skips every file that can't contain the key. Akamai's CDN used one to cache an object only on its second request, because most objects are requested exactly once, as Bruce Maggs and Ramesh Sitaraman describe in "Algorithmic nuggets in content delivery". In both cases a false positive costs one wasted lookup, and a false negative, which never happens, would have cost a wrong answer.

Count-min sketch: frequencies that only overcount

A count-min sketch, from Graham Cormode and S. Muthukrishnan in 2005, is d rows of w counters, each row with its own hash function. Adding an item increments one counter in each row. The estimate for an item is the smallest of its d counters. Every counter the item touches holds its true count plus whatever else hashed to the same place, so each one is an overestimate, and the minimum is the least wrong. With w = ⌈e/ε⌉ and d = ⌈ln(1/δ)⌉, the estimate exceeds the true count by at most εN, where N is the total number of events, with probability at least 1 − δ.

That bound is relative to the whole stream, not to the item's own count. An error of εN is negligible for an item seen thousands of times and enormous for one seen twice. So the count-min sketch is a tool for heavy hitters: the most frequent source addresses on a network link, trending search terms, the noisiest API keys for rate limiting. Combined with a small heap of the current leaders (chapter 3), it answers "what are the top 100?" over a stream it never stores.

HyperLogLog: counting distinct items in kilobytes

Start with coins. If someone tells you the longest run of heads they saw before a tail was 10, they probably flipped about 2¹⁰ times, since a run that long is a one-in-a-thousand event. Hash every item to a random-looking 64-bit number and look at its leading zero bits: a hash starting with r − 1 zeros and then a one turns up with probability 2^−r, so the largest such r seen estimates log₂ of the number of distinct items. Duplicates hash identically and can't raise the maximum, which is exactly why the estimate counts distinct items.

A single maximum is noisy, since one lucky hash could claim a billion items. HyperLogLog, from Philippe Flajolet and colleagues in 2007, splits the hashes among m registers using their first p bits, keeps the largest rank seen in each register, and combines the registers with a harmonic mean, which damps the occasional huge value. Its standard error is 1.04/√m: 3.25% with 1,024 registers, 0.81% with 16,384. A rank never exceeds 64, so a register needs only 6 bits, and 16,384 registers fit in 12 KB, which is what Redis's PFCOUNT uses. The example spends a whole byte per register for clarity. For small counts, when many registers are still empty, it switches to linear counting, an estimate from the number of empty registers that's more accurate in that range.

All three, measured

The example builds all three from the definitions above and measures them against exact answers: a Bloom filter at two sizes, a count-min sketch over a skewed stream in which a few items are very common, and HyperLogLog at two register counts.

import hashlib, math, random

def h128(item):
    d = hashlib.blake2b(str(item).encode(), digest_size=16).digest()
    return int.from_bytes(d[:8], "little"), int.from_bytes(d[8:], "little") | 1

class BloomFilter:
    """m bits and k hash functions: 'definitely not present' or 'probably present'."""
    def __init__(self, m, k):
        self.m, self.k, self.bits = m, k, bytearray((m + 7) // 8)
    def _positions(self, item):
        a, b = h128(item)
        return [(a + i * b) % self.m for i in range(self.k)]   # k hashes from two (Kirsch-Mitzenmacher)
    def add(self, item):
        for p in self._positions(item): self.bits[p // 8] |= 1 << (p % 8)
    def __contains__(self, item):
        return all(self.bits[p // 8] >> (p % 8) & 1 for p in self._positions(item))

class CountMin:
    """d rows of w counters; an item's estimate is its smallest counter, never an undercount."""
    def __init__(self, eps, delta):
        self.w, self.d = math.ceil(math.e / eps), math.ceil(math.log(1 / delta))
        self.rows = [[0] * self.w for _ in range(self.d)]
    def _cols(self, item):
        a, b = h128(item)
        return [(a + r * b) % self.w for r in range(self.d)]
    def add(self, item):
        for r, c in enumerate(self._cols(item)): self.rows[r][c] += 1
    def estimate(self, item):
        return min(self.rows[r][c] for r, c in enumerate(self._cols(item)))

class HyperLogLog:
    """Estimates how many distinct items it has seen from the longest runs of zero bits in their hashes."""
    def __init__(self, p):
        self.p, self.m, self.reg = p, 1 << p, bytearray(1 << p)
    def add_hash(self, x):
        j = x >> (64 - self.p)                                  # the first p bits pick a register
        rest = (x << self.p) & ((1 << 64) - 1)
        rank = 64 - self.p + 1 if rest == 0 else (64 - rest.bit_length()) + 1   # leading zeros + 1
        self.reg[j] = max(self.reg[j], rank)
    def count(self):
        alpha = 0.7213 / (1 + 1.079 / self.m)
        raw = alpha * self.m * self.m / sum(2.0 ** -r for r in self.reg)
        zeros = self.reg.count(0)
        if raw <= 2.5 * self.m and zeros:                      # small-range correction: linear counting
            return self.m * math.log(self.m / zeros)
        return raw

random.seed(7)
members = [f"user-{i}" for i in range(10_000)]
others = [f"visitor-{i}" for i in range(100_000)]
print("Bloom filter, 10,000 members, 100,000 non-members tested:")
for bits_per_key in (8, 12):
    m = bits_per_key * len(members)
    k = round(bits_per_key * math.log(2))                       # the k that minimises false positives
    bf = BloomFilter(m, k)
    for x in members: bf.add(x)
    assert all(x in bf for x in members)                        # no false negatives, ever
    fp = sum(x in bf for x in others) / len(others)
    theory = (1 - math.exp(-k * len(members) / m)) ** k
    print(f"   {bits_per_key:2} bits per key, k={k}: false positives {fp:.2%} (formula {theory:.2%}), "
          f"{m // 8:,} bytes")

ranks = list(range(1, 10_001))
weights = [1 / r for r in ranks]                                # a few items are very common
stream = random.choices(ranks, weights, k=200_000)
truth = {}
for x in stream: truth[x] = truth.get(x, 0) + 1
cm = CountMin(eps=0.001, delta=0.01)
for x in stream: cm.add(x)
errors = [cm.estimate(x) - truth[x] for x in truth]
print(f"count-min sketch ({cm.d} x {cm.w} = {cm.d * cm.w:,} counters) over {len(stream):,} events, {len(truth):,} items:")
print(f"   item 1: true {truth[1]:,}, estimate {cm.estimate(1):,}   item 500: true {truth[500]}, estimate {cm.estimate(500)}")
print(f"   never undercounts: {min(errors) >= 0};  worst overcount {max(errors)} (bound eps*N = {0.001 * len(stream):.0f}); "
      f"{sum(e <= 0.001 * len(stream) for e in errors) / len(errors):.1%} of items within it")

stream = [random.randrange(300_000) for _ in range(1_000_000)]
true_distinct = len(set(stream))
sketches = {p: HyperLogLog(p) for p in (10, 14)}
for x in stream:
    a, _ = h128(x)
    for s in sketches.values(): s.add_hash(a)
print(f"HyperLogLog over {len(stream):,} events with {true_distinct:,} distinct values:")
for p, s in sketches.items():
    est = s.count()
    print(f"   {s.m:6,} registers ({s.m:,} bytes): estimate {est:,.0f}, error {abs(est - true_distinct) / true_distinct:.2%} "
          f"(typical error 1.04/sqrt(m) = {1.04 / math.sqrt(s.m):.2%})")

Bloom filter, 10,000 members, 100,000 non-members tested:
    8 bits per key, k=6: false positives 2.14% (formula 2.16%), 10,000 bytes
   12 bits per key, k=8: false positives 0.28% (formula 0.31%), 15,000 bytes
count-min sketch (5 x 2719 = 13,595 counters) over 200,000 events, 9,653 items:
   item 1: true 20,345, estimate 20,350   item 500: true 38, estimate 46
   never undercounts: True;  worst overcount 96 (bound eps*N = 200); 100.0% of items within it
HyperLogLog over 1,000,000 events with 289,162 distinct values:
    1,024 registers (1,024 bytes): estimate 283,435, error 1.98% (typical error 1.04/sqrt(m) = 3.25%)
   16,384 registers (16,384 bytes): estimate 289,029, error 0.05% (typical error 1.04/sqrt(m) = 0.81%)

The Bloom filter matches the formula closely. At 8 bits per key with k = 6 it measured 2.14% false positives against 2.16% predicted; at 12 bits per key with k = 8, 0.28% against 0.31%. Every one of the 10,000 members was found. The whole filter is 10,000 bytes, one byte per key, for keys that are 6 to 9 characters long before any hash table overhead.

The count-min sketch shows the shape of its error. The most common item, seen 20,345 times, was estimated at 20,350, off by 0.02%. Item 500, seen 38 times, was estimated at 46, off by 21%. The worst overcount among all 9,653 items was 96, inside the bound of 200, and nothing was undercounted. It does this in 13,595 counters, whatever the length of the stream.

HyperLogLog estimated 289,162 distinct values to within 1.98% using 1,024 registers and 0.05% using 16,384. The second run was luckier than its typical 0.81%, because the standard error describes a spread, not a limit. An exact count would have to store all 289,162 values, at least 2.3 MB at 8 bytes each before any table overhead, and that grows with the data. The sketch stays at 16 KB, or 12 KB packed, at a billion distinct values.

Choosing, sizing and merging

Merging matters as much as memory. Two sketches with the same size and hash functions, built on different machines or on different days, combine into exactly the sketch of the combined data. That's how analytics systems count distinct values across shards without moving the values: BigQuery's APPROX_COUNT_DISTINCT and Trino's approx_distinct are HyperLogLog. Exact sets or counts would have to be shipped in full.

Size for the peak, not the average. The false-positive formula assumed you knew n. A filter sized for one number of keys and then given twice as many is no longer the filter you designed (exercise 7.1). A sketch sized for today's traffic needs headroom for next year's.

Easy to check, hard to find

Take subset sum: given a list of numbers and a target, is there a subset that adds up to exactly the target? If someone hands you a subset, checking it takes a few additions. Finding one is another matter. The obvious method tries all 2ⁿ subsets, and nobody knows a method that is polynomial in the size of the input for every instance.

That asymmetry defines the class. P is the set of problems that can be solved in polynomial time. NP is the set of yes-or-no problems whose yes-answers come with a certificate, like the subset, that can be checked in polynomial time. The name means "nondeterministic polynomial": a machine that could guess the right certificate would finish in polynomial time. Every problem in P is also in NP, since you can check an answer by solving the problem yourself. Whether P equals NP is the most famous open question in computer science. Stephen Cook posed it in 1971, it's one of the Clay Millennium Prize problems, and most researchers believe the answer is no.

Reductions: why one problem stands for thousands

A reduction from problem A to problem B is a polynomial-time translation of any instance of A into an instance of B with the same answer. Given one, a fast algorithm for B would give a fast algorithm for A, so B is at least as hard as A. Cook, and independently Leonid Levin, proved that Boolean satisfiability (SAT) is as hard as everything in NP: any NP problem can be translated into a formula that's satisfiable exactly when the answer is yes. In 1972 Richard Karp reduced SAT to 21 classic problems, among them vertex cover, clique, Hamiltonian cycle, graph colouring and a version of subset sum. Garey and Johnson's 1979 book catalogued hundreds more, and there are now thousands. A problem is NP-hard if everything in NP reduces to it, and NP-complete if it's also in NP itself.

Direction is the classic mistake. To show that your new problem X is hard, reduce a known NP-complete problem to X: "if I could solve X quickly, I could solve vertex cover quickly." Reducing X to SAT shows only that X is no harder than SAT. That's useful too, in a different way, because it means a SAT solver can solve X.

Feeling the wall

The example makes the wall concrete. It brute-forces subset sum on instances with no solution, which is the worst case because every subset has to be tried. It checks a certificate, solves a larger instance by dynamic programming, solves a 10-city travelling salesperson problem exactly and with two heuristics, and compares an approximation for vertex cover against the true optimum.

import itertools, math, random

random.seed(8)
# 1. Subset sum: checking an answer is fast; finding one, in general, means trying subsets.
def brute_force(nums, target):
    """Sum of every subset, built incrementally: 2^n subsets, each one addition."""
    sums = [0] * (1 << len(nums))
    for mask in range(1, 1 << len(nums)):
        low = mask & -mask
        sums[mask] = sums[mask ^ low] + nums[low.bit_length() - 1]
        if sums[mask] == target:
            return mask, mask
    return None, len(sums) - 1

print("subset sum, no subset hits the target (the worst case), subsets checked:")
for n in (10, 14, 18, 20):
    nums = [random.randrange(1, 10**9) * 2 for _ in range(n)]  # even numbers can never sum to an odd target
    found, checked = brute_force(nums, 1)
    print(f"   n = {n:2}: {checked:>17,} subsets examined, answer found: {found is not None}")
print(f"   n = 40: {2**40 - 1:>17,} subsets: each extra number doubles the work")

nums = [random.randrange(1, 1000) for _ in range(30)]
certificate = random.sample(range(30), 12)
target = sum(nums[i] for i in certificate)
print(f"verifying a claimed answer for n = 30 takes {len(certificate)} additions: "
      f"{sum(nums[i] for i in certificate) == target}  (brute force could need 2^30 = {2**30:,} subsets)")

reachable, work = {0}, 0
for x in nums:                                                  # dynamic programming over sums, not subsets
    work += len(reachable)
    reachable |= {s + x for s in reachable if s + x <= target}
print(f"dynamic programming with small numbers (all below 1,000): target {target:,} reachable: "
      f"{target in reachable}, {work:,} steps (bound: 30 x {target + 1:,} sums = {30 * (target + 1):,})")
print(f"the same idea with 40-bit numbers would need tables of about 2^40 = {2**40:,} cells")

# 2. Travelling salesperson: exact by brute force, then two fast heuristics.
pts = [(random.random(), random.random()) for _ in range(10)]
def length(tour): return sum(math.dist(pts[tour[i]], pts[tour[i - 1]]) for i in range(len(tour)))
best = min(((0,) + p for p in itertools.permutations(range(1, 10))), key=length)
nn, left = [0], set(range(1, 10))
while left:                                                     # nearest neighbour: always go to the closest city
    nxt = min(left, key=lambda c: math.dist(pts[nn[-1]], pts[c])); nn.append(nxt); left.remove(nxt)
tour, improved = nn[:], True
while improved:                                                 # 2-opt: reverse any segment that shortens the tour
    improved = False
    for i in range(1, len(tour) - 1):
        for j in range(i + 1, len(tour)):
            new = tour[:i] + tour[i:j + 1][::-1] + tour[j + 1:]
            if length(new) < length(tour) - 1e-12:
                tour, improved = new, True
print(f"\ntravelling salesperson, 10 cities ({math.factorial(9):,} tours tried for the exact answer):")
print(f"   optimal {length(best):.3f}   nearest neighbour {length(nn):.3f} (+{length(nn) / length(best) - 1:.0%})   "
      f"after 2-opt {length(tour):.3f} (+{length(tour) / length(best) - 1:.0%})")

# 3. Vertex cover: an approximation with a guarantee.
nodes = range(14)
edges = random.sample([(a, b) for a in nodes for b in nodes if a < b], 25)
def covers(s): return all(a in s or b in s for a, b in edges)
exact = min((set(c) for r in range(15) for c in itertools.combinations(nodes, r) if covers(set(c))), key=len)
approx = set()
for a, b in edges:                                              # take both ends of any edge still uncovered
    if a not in approx and b not in approx:
        approx |= {a, b}
print(f"\nvertex cover, 14 nodes and 25 edges: optimal {len(exact)} nodes (found by checking all 2^14 = {2**14:,} subsets)")
print(f"   matching-based approximation: {len(approx)} nodes, covers every edge: {covers(approx)}, "
      f"within 2x of optimal: {len(approx) <= 2 * len(exact)}")

subset sum, no subset hits the target (the worst case), subsets checked:
   n = 10:             1,023 subsets examined, answer found: False
   n = 14:            16,383 subsets examined, answer found: False
   n = 18:           262,143 subsets examined, answer found: False
   n = 20:         1,048,575 subsets examined, answer found: False
   n = 40: 1,099,511,627,775 subsets: each extra number doubles the work
verifying a claimed answer for n = 30 takes 12 additions: True  (brute force could need 2^30 = 1,073,741,824 subsets)
dynamic programming with small numbers (all below 1,000): target 6,575 reachable: True, 102,145 steps (bound: 30 x 6,576 sums = 197,280)
the same idea with 40-bit numbers would need tables of about 2^40 = 1,099,511,627,776 cells

travelling salesperson, 10 cities (362,880 tours tried for the exact answer):
   optimal 3.071   nearest neighbour 3.201 (+4%)   after 2-opt 3.108 (+1%)

vertex cover, 14 nodes and 25 edges: optimal 8 nodes (found by checking all 2^14 = 16,384 subsets)
   matching-based approximation: 10 nodes, covers every edge: True, within 2x of optimal: True

Each extra number doubles the brute-force work. Twenty numbers mean a million subsets, which is nothing. Forty mean a trillion, which compiled code can still grind through. Sixty mean a million times more again. Checking a claimed answer for 30 numbers took 12 additions, against a search space of more than a billion subsets. That gap between finding and checking is NP in one line of output.

What to do instead

Check the size first. Exponential doesn't automatically mean too slow. If your instances have 20 items, brute force over a million subsets is the simplest correct solution. With pruning (backtracking, branch and bound), brute force often reaches much further in practice than its worst case suggests.

Exploit small numbers. The dynamic program tracks which sums are reachable instead of which subsets exist. For 30 numbers below 1,000 it decided reachability of 6,575 in 102,145 steps. That's polynomial in the value of the target, not in its length in bits, which makes it pseudo-polynomial. With 40-bit numbers the table would need about 2⁴⁰ cells and the trick is gone. Problems like subset sum and knapsack are called weakly NP-complete for this reason: they're practical whenever the numbers involved are small.

Exploit structure. NP-hardness is a statement about the worst instances, and yours may avoid them. Vertex cover is easy on trees and on bipartite graphs, where by Kőnig's theorem it equals a maximum matching. 2-SAT is solvable in linear time though 3-SAT is NP-complete. Colouring an interval graph, such as assigning meetings to rooms, is solved by a greedy algorithm. Some problems are hard only in one parameter: a vertex cover of size k can be found in about 2^k·n steps by branching on the two ends of an uncovered edge, which is fast for a huge graph with a small cover.

Approximate with a guarantee. The vertex-cover algorithm in the example takes both ends of any uncovered edge. It found 10 nodes where the optimum is 8, and it can never be worse than twice the optimum (exercise 8.3). For travelling salesperson with real distances, which obey the triangle inequality, Christofides' algorithm guarantees a tour within 1.5 times the optimum, an approximation ratio improved in 2020 by a famously tiny margin. Some problems resist even this: without the triangle inequality, no polynomial algorithm can approximate travelling salesperson within any constant factor unless P = NP.

Use heuristics, and measure them. Nearest neighbour built a tour 4% longer than the optimum; 2-opt, which reverses any segment that shortens the tour, brought it to 1%. Neither guarantees anything in general, but local search of this kind routinely gets within a few percent on real instances. Measure a heuristic against exact answers on small instances, or against a lower bound on large ones, before you trust it.

Hand it to a solver. SAT solvers and mixed-integer programming solvers routinely solve industrial instances with millions of variables, in hardware verification, scheduling and logistics. Their worst case is still exponential, but real instances have structure that decades of solver engineering exploit. Translating your problem into SAT or an integer program, which is itself a reduction, and running Google's OR-Tools, the open-source HiGHS or a commercial solver like Gurobi often beats a hand-written algorithm by a wide margin.

Hard problems in disguise

NP-hard problems rarely announce themselves. Choosing one version of each package so that every dependency constraint holds is NP-complete: Russ Cox showed how to encode 3-SAT in it. That's why Go modules chose minimal version selection, which sidesteps the general problem by design, while Dart's pub and Python's uv use PubGrub, a resolver built on SAT-solving techniques. Scheduling exams so no student has a clash is graph colouring, and so is register allocation in a compiler. Placing containers onto servers is bin packing. Choosing which features fit a release budget is knapsack. Planning a delivery route is travelling salesperson. When a problem asks you to choose a subset, an assignment or an ordering that satisfies constraints or optimises a total, suspect NP-hardness before you spend a week looking for a clever polynomial algorithm.