Randomized Algorithms
Probability as a Tool

History: Monte Carlo method (Ulam, von Neumann, Metropolis, 1946–49). Rabin (1976) and Solovay–Strassen (1977) primality tests started the modern field.

n = 1_000_000
arr = [1] * (n // 2) + [0] * (n // 2); rng.shuffle(arr)

def lv_find(a):                  # Las Vegas: always right, time random
    probes = 0
    while True:
        i = rng.randrange(len(a)); probes += 1
        if a[i] == 1: return i, probes

def mc_find(a, k):               # Monte Carlo: time fixed, may fail
    for _ in range(k):
        i = rng.randrange(len(a))
        if a[i] == 1: return i
    return None

pr = [lv_find(arr)[1] for _ in range(10_000)]
fails = sum(mc_find(arr, 10) is None for _ in range(100_000))

LV: mean probes 2.008, worst 14
MC k=10: fail rate 0.00100, bound 2^-10 = 0.00098

# 1000 fair coins, 20000 trials: how often |heads - 500| >= 50?
dev = sum(abs(sum(rng.random() < .5 for _ in range(1000)) - 500) >= 50
          for _ in range(20000)) / 20000
cheb  = (1000 * .25) / 50**2              # Var = n/4
chern = 2 * math.exp(-2 * 50**2 / 1000)   # Hoeffding

Pr[|X-500|>=50]: measured 0.0018, Chebyshev <= 0.100, Hoeffding <= 0.0135

Step 1. Let X_i = 1 if run i is wrong. So E[X] = k/3, where X = ΣX_i.

Step 2. The majority is wrong only if X ≥ k/2. That is t = k/2 − k/3 = k/6 above the mean.

Step 3. Apply Hoeffding (one-sided):

Pr[X ≥ k/2] ≤ e^−2t²/k = e^{−2(k/6)²/k} = e^−k/18

Step 4. Solve e^−k/18 ≤ 10⁻⁶:

k ≥ 18 · ln(10⁶) = 248.7 ⇒ k = 249

p = 2/3
def majority_error(k):           # exact binomial tail, k odd
    return sum(comb(k, i) * p**i * (1-p)**(k-i) for i in range(k//2 + 1))
print("  k   exact     bound e^(-k/18)")
for k in (1, 9, 49, 99, 193, 249):
    print(f"{k:>3}   {majority_error(k):.2e}  {math.exp(-k/18):.2e}")

  k   exact     bound e^(-k/18)
  1   3.33e-01  9.46e-01
  9   1.45e-01  6.07e-01
 49   7.87e-03  6.57e-02
 99   3.09e-04  4.09e-03
193   9.03e-07  2.20e-05
249   2.96e-08  9.82e-07

Common mistake: reusing the same seed on each run. Then the runs are not independent, and Chernoff does not apply. Also pick an odd k, or set a tie rule.

Proof. Name the keys by sorted rank: z₁ < z₂ < … < z_n. Let X_ij = 1 if z_i and z_j are ever compared.

• Two keys are compared only if one is the pivot. After that the pivot is gone. So each pair is compared at most once: C = Σ_i<j X_ij.
• Look at the block Z_ij = {z_i, …, z_j}. It stays together until the first pivot from inside it is chosen.
• If that first pivot is z_i or z_j, they get compared. If it is anything between, they are split forever.
• Each of the j−i+1 keys is equally likely to be first: Pr[X_ij=1] = 2/(j−i+1).

E[C] = Σ_i<j 2/(j−i+1) = Σ_d=2ⁿ (n−d+1)·2/d

= 2(n+1)H_n − 4n  □

Sort [7, 3, 9, 1, 5, 8, 2]. Each pivot is picked at random. Here is one run. Every non-pivot is compared once with the pivot.

call [7, 3, 9, 1, 5, 8, 2]   pivot 5   6 compares
     less [3, 1, 2]    more [7, 9, 8]
call [3, 1, 2]               pivot 1   2 compares
     less []           more [3, 2]
call [3, 2]                  pivot 3   1 compare
     less [2]          more []
call [7, 9, 8]               pivot 8   2 compares
     less [7]          more [9]
                                total 11 compares

Try it: redo the trace with pivots 2, then 8, then 5. (Answer: 6 + 4 + 2 = 12 compares.)

def rquicksort(a, counter):
    if len(a) <= 1: return a
    p = a[rng.randrange(len(a))]
    counter[0] += len(a) - 1      # pivot vs every other item
    less = [x for x in a if x < p]
    eq   = [x for x in a if x == p]
    more = [x for x in a if x > p]
    return rquicksort(less, counter) + eq + rquicksort(more, counter)

def H(n): return sum(1 / k for k in range(1, n + 1))
# mean of 20 runs per n, input already sorted: list(range(n))

n      measured    2(n+1)H_n-4n
100           629            648
1000        11131          10986
10000      154135         155772
sorted input n=2000 comparisons: 25328

Sorted input is the classic worst case for "pivot = first element": n(n−1)/2 = 1,999,000 comparisons at n=2000. The random pivot needs 25,328, about 79× fewer.

def quickselect(a, k):
    """k-th smallest (0-based). Expected O(n)."""
    while True:
        p = a[rng.randrange(len(a))]
        less = [x for x in a if x < p]
        eq   = [x for x in a if x == p]
        if k < len(less): a = less
        elif k < len(less) + len(eq): return p
        else:
            k -= len(less) + len(eq)
            a = [x for x in a if x > p]

data = [rng.randint(0, 10**6) for _ in range(10001)]
print("median:", quickselect(data, 5000) == sorted(data)[5000])

median: True

def karger(edges, n):
    parent = list(range(n))
    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]; x = parent[x]
        return x
    order = edges[:]; rng.shuffle(order)   # random order = random contractions
    comps = n
    for u, v in order:
        if comps == 2: break
        ru, rv = find(u), find(v)
        if ru != rv: parent[ru] = rv; comps -= 1
    return sum(find(u) != find(v) for u, v in edges)

Graph from the last slide: edges ab, ac, bc, bd, cd. Vertices a and d have degree 2. So the min cut is 2, and there are two of them: {a}|{bcd} and {d}|{abc}.

Common mistake: deleting parallel edges after a merge. Keep them. They make a heavy link more likely to be contracted, which is what protects a small cut.

Proof. Let k = |C|. Every vertex has degree ≥ k, or it would be a smaller cut by itself. So the graph has m ≥ nk/2 edges.

At step i, with r = n − i + 1 super-vertices left, the same holds: m_i ≥ rk/2. So:

Pr[pick an edge of C] ≤ k/(rk/2) = 2/r

Pr[C survives] ≥ ∏_r=3ⁿ (1 − 2/r) = ∏ (r−2)/r

= (1/3)(2/4)(3/5)…((n−2)/n) = 2/(n(n−1))  □

The product telescopes. A bonus corollary: any graph has at most C(n,2) distinct min cuts.

# two 5-cliques joined by 2 edges: min cut = 2
E, n = two_cliques(), 10
hits = sum(karger(E, n) == 2 for _ in range(10000))
T = math.ceil(n * (n - 1) / 2 * math.log(100))
best = min(karger(E, n) for _ in range(T))

single run success: 0.274   bound 2/(n(n-1)) = 0.022
repeat 208 times -> min cut 2

Karger–Stein (1996): contract only down to n/√2, then recurse twice. Success is Ω(1/log n), total time O(n² log³ n).

Next: algorithms that check instead of compute. Primes, matrices, strings, polynomials.

def is_probable_prime(n, k=20):
    if n < 2: return False
    for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29):
        if n % p == 0: return n == p
    d, s = n - 1, 0
    while d % 2 == 0: d //= 2; s += 1   # n-1 = 2^s * d, d odd
    for _ in range(k):
        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 composite
    return True

print(is_probable_prime(2**127 - 1), is_probable_prime(561),
      is_probable_prime(2**127 + 1))

True False False

Tested against a sieve for every n < 100,000: no mistakes.

n=561: strong liars 8/558 = 0.014
n=1105: strong liars 28/1102 = 0.025
n=8911: strong liars 1780/8908 = 0.200
n=221: strong liars 4/218 = 0.018
Fermat liars for 561: 318 of 558

561 = 3·11·17 is a Carmichael number. It fools the plain Fermat test for every base coprime to it, but Miller–Rabin catches it. Error is one-sided: "composite" is always right.

2^35  mod 561 = 263
2^70  mod 561 = 166     (263^2 mod 561)
2^140 mod 561 = 67      (166^2 mod 561)
2^280 mod 561 = 1       (67^2  mod 561)
2^560 mod 561 = 1       Fermat test: "looks prime"

Bonus: gcd(67−1, 561) = 33 and gcd(67+1, 561) = 17. The witness even hands us factors: 561 = 3 · 11 · 17.

5^3  mod 97 = 28
5^6  mod 97 = 8
5^12 mod 97 = 64
5^24 mod 97 = 22
5^48 mod 97 = 96   = n-1, i.e. -1  -> pass

Common mistake: computing 2⁵⁶⁰ first, then reducing. Use pow(a, d, n). It reduces at every step.

def matvec(A, x):
    return [sum(a * b for a, b in zip(row, x)) for row in A]

def freivalds(A, B, C, k=20):
    n = len(A)
    for _ in range(k):
        r = [rng.randint(0, 1) for _ in range(n)]
        if matvec(A, matvec(B, r)) != matvec(C, r):
            return False          # definitely wrong
    return True                   # right with prob >= 1 - 2^-k

n=120: multiply 0.115s, Freivalds k=20 0.043s, correct C -> True
one-entry error caught by a single round: 0.514

Each round costs 3 matrix-vector products, O(n²). The naive triple-loop multiply costs O(n³), and the gap widens as n grows.

Rolling hash: treat a window as a number in base b, taken mod q:

h(s_i+1..i+m) = (h(s_i..i+m−1) − s_ib^m−1)·b + s_i+m

Sliding the window costs O(1). Expected total time is O(n + m).

def rabin_karp(text, pat, base=256, mod=(1 << 61) - 1):
    m, n = len(pat), len(text)
    if m > n: return []
    hp = ht = 0
    for i in range(m):
        hp = (hp * base + ord(pat[i])) % mod
        ht = (ht * base + ord(text[i])) % mod
    top = pow(base, m - 1, mod)
    hits = []
    for i in range(n - m + 1):
        if hp == ht and text[i:i + m] == pat:  # verify: Las Vegas
            hits.append(i)
        if i + m < n:
            ht = ((ht - ord(text[i]) * top) * base
                  + ord(text[i + m])) % mod
    return hits

print(rabin_karp("abracadabra" * 3, "abra"))

[0, 7, 11, 18, 22, 29]

With the verify step it is Las Vegas: never wrong, only slower on hash collisions. Drop the verify and it becomes Monte Carlo. Also checked against str.startswith on 200 random strings.

The same rolling-hash trick powers rsync block matching and content-defined chunking in backup tools.

P = (1 << 61) - 1          # a Mersenne prime, field Z_P
def poly_equal(f, g, nvars, trials=10):
    for _ in range(trials):
        xs = [rng.randrange(P) for _ in range(nvars)]
        if f(*xs) % P != g(*xs) % P:
            return False
    return True

f = lambda x, y: (x + y) ** 5
g = lambda x, y: (x**5 + 5*x**4*y + 10*x**3*y**2
                  + 10*x**2*y**3 + 5*x*y**4 + y**5)
h = lambda x, y: g(x, y) + x*y            # off by one term
print(poly_equal(f, g, 2), poly_equal(f, h, 2))

True False

Error per trial is ≤ 5/2⁶¹. The test file also checks the 4×4 Vandermonde determinant identity this way.

p_ = 2**31 - 1
def make_hash(m):
    a, b = rng.randrange(1, p_), rng.randrange(p_)
    return lambda x: ((a * x + b) % p_) % m

Pr[h(x)=h(y)] measured 0.0100  vs 1/m = 0.0100

class SkipList:
    MAXH = 32
    def __init__(self):
        self.head = Node(None, self.MAXH); self.h = 1
    def _height(self):                  # Pr[h >= i] = 2^-(i-1)
        h = 1
        while h < self.MAXH and rng.random() < 0.5: h += 1
        return h
    def insert(self, key):
        update, x = [self.head] * self.MAXH, self.head
        for lvl in reversed(range(self.h)):       # right, then down
            while x.next[lvl] and x.next[lvl].key < key: x = x.next[lvl]
            update[lvl] = x
        h = self._height(); self.h = max(self.h, h); node = Node(key, h)
        for lvl in range(h):                      # splice into h lists
            node.next[lvl] = update[lvl].next[lvl]; update[lvl].next[lvl] = node
    # search: same walk as insert, counting steps right and down

N=  1000: height 10, mean search steps 19.0, 2*log2(N) = 19.9
N=100000: height 18, mean search steps 32.2, 2*log2(N) = 33.2

Why 2 log₂ n? Walk the path backwards. At each node, go up w.p. 1/2 or left w.p. 1/2. So about 2 steps per level, and about log₂ n levels.

Used in Redis sorted sets, LevelDB/RocksDB memtables, and Java's ConcurrentSkipListMap. Easy to make lock-free.

class Bloom:
    def __init__(self, m, k):
        self.m, self.k, self.bits = m, k, bytearray(m)
    def _idx(self, item):
        h = hashlib.sha256(item.encode()).digest()
        h1 = int.from_bytes(h[:8], "big")
        h2 = int.from_bytes(h[8:16], "big") | 1
        return [(h1 + i * h2) % self.m for i in range(self.k)]
    def add(self, item):
        for i in self._idx(item): self.bits[i] = 1
    def __contains__(self, item):
        return all(self.bits[i] for i in self._idx(item))

nitems, mbits = 10000, 100000
k = round(mbits / nitems * math.log(2))           # 7
bf = Bloom(mbits, k)
for i in range(nitems): bf.add(f"user{i}")
fp = sum(f"other{i}" in bf for i in range(100000)) / 100000

m/n = 10, k = 7: measured FP 0.0084, formula 0.0082

Double hashing h₁ + i·h₂ gives k indices from one digest (Kirsch–Mitzenmacher 2006). All 10,000 inserted keys test positive: no false negatives.

Uses: LSM-tree lookups (RocksDB, Cassandra), CDN "one-hit wonder" filters, Chrome's old Safe Browsing list. Variants: counting Bloom filters (allow deletes), cuckoo and xor filters.

Problem. A stream of unknown length goes by once. Keep a uniform random sample of k items using O(k) memory.

def reservoir(stream, k):          # Algorithm R (Vitter 1985)
    R = []
    for i, x in enumerate(stream):
        if i < k:
            R.append(x)
        else:
            j = rng.randrange(i + 1)   # 0..i
            if j < k:
                R[j] = x               # keep x with prob k/(i+1)
    return R

counts = [0] * 10
for _ in range(50000):
    for x in reservoir(range(10), 3): counts[x] += 1
print(counts)                          # expect 15000 each

[14818, 14899, 15001, 15095, 15047, 14681, 15058, 15104, 15167, 15130]

Each count should be 50000 × 3/10 = 15000. The standard deviation is about 102. Nine counts sit within 1.5σ. One (14681) is 3.1σ low: unlucky, not a bug. With ten counts, one outlier like that turns up now and then. Rerun with other seeds before you suspect the code.

for N in (10**3, 10**4, 10**5, 10**6):
    inside = sum(rng.random()**2 + rng.random()**2 < 1
                 for _ in range(N))
    est = 4 * inside / N
    print(f"{N:<9} {est:.5f}   {abs(est - math.pi):.5f}")

N         estimate    error
1000      3.21600   0.07441
10000     3.14800   0.00641
100000    3.13116   0.01043
1000000   3.14318   0.00159

Each point lands in the quarter circle with probability π/4. The estimator 4·inside/N is unbiased.

Error is not monotone: N = 10⁵ was worse than 10⁴ in this run. Randomness only promises the typical error.

def coupons(n):
    seen, t = set(), 0
    while len(seen) < n:
        seen.add(rng.randrange(n)); t += 1
    return t

coupons n=10: mean 29.2, n*H_n = 29.3
coupons n=100: mean 518.5, n*H_n = 518.7
coupons n=1000: mean 7497.1, n*H_n = 7485.5

Shows up in load balancing: n balls into n bins leaves about n/e bins empty. Filling every bin needs n ln n balls.

def birthday_trial(n_people, days=365):
    return len({rng.randrange(days) for _ in range(n_people)}) < n_people

10 people: simulated 0.117, exact 0.117
23 people: simulated 0.501, exact 0.507
50 people: simulated 0.970, exact 0.970

def walk_2sat(cnf, n, max_steps=None):
    max_steps = max_steps or 100 * n * n
    a = {v: rng.random() < .5 for v in range(1, n + 1)}
    for step in range(max_steps):
        bad = [c for c in cnf
               if not any(a[abs(l)] == (l > 0) for l in c)]
        if not bad: return a, step
        l = rng.choice(rng.choice(bad))   # random literal, random bad clause
        a[abs(l)] = not a[abs(l)]
    return None, max_steps

# 50 random satisfiable 2-CNFs, n=50 vars, 150 clauses
# (clauses sampled to agree with a hidden assignment)

2-SAT walk n=50: mean flips 141, max 416, budget 2n^2 = 5000

Each flip here rescans all clauses, which is fine for a demo. A real version keeps a list of unsatisfied clauses.

Schöning (1999) used the same idea for 3-SAT. Restart every 3n flips and it runs in O((4/3)ⁿ), far better than 2ⁿ. WalkSAT is the practical heir.