46Randomization & the probabilistic mind

In Chapter 45 we bent strings and number theory into the machinery of crypto and quant. And every algorithm we've built, all volume long, has made its choices deterministically: the same pivot, the same hash, the same next node, every single time. In this chapter we do the one thing that sounds like surrender and turns out to be a superpower. We flip a coin. Here's the plan. We'll armor a sort so no attacker can find its worst case. We'll compute a number like π by throwing darts at a wall. We'll price a financial option by simulating a million possible futures. And we'll draw a perfectly fair sample from a stream far too big to hold. The whole way through we keep asking the one question that actually matters: how can deliberately adding randomness make an algorithm more reliable, not less? By the end you'll reach for a coin on reflex. You'll reach for it when being predictable is the weakness, when the exact answer is out of reach, and when you need fairness over data you can't even count.

01Every fixed plan has a nemesis

Let's start with a word we've leaned on all volume without ever saying it out loud: deterministic. A deterministic algorithm always makes the same choice on the same input. Feed it the same list twice and it does the exact same thing, step for step. Usually that's a virtue, because it's what makes code testable and predictable. But watch what it costs us. A deterministic algorithm has a fixed worst case — one specific arrangement of input that drives it to its slowest possible behavior, every time. And suppose someone can read your code, because it's open-source, or a public API, or a well-known library. Then they can construct that exact arrangement and hand it to you on purpose.

Quicksort is the classic victim of exactly this attack. Tell it "always pick the first element as the pivot," then hand it an already-sorted list, and watch it unravel. Every partition peels off just one element, so the recursion sinks n levels deep. That drops the tidy O(n log n) straight into O(n²). You stepped quicksort's partition in trace T28, and here we attack that same code. This is not a thought experiment, either — it's a working denial-of-service. Put a naive quicksort behind a web endpoint, and an attacker who sends sorted data can pin your CPU flat with a tiny request.

The escape is not a cleverer fixed rule, because the adversary just studies the new rule and defeats that one too. The escape is to stop being predictable altogether. Once the algorithm flips a coin to make its choices, no single input is reliably its worst case. The worst case now depends on coins the attacker cannot see.

Before we lean on coin flips, let's be honest about where they come from. Python's random module doesn't flip physical coins, because it runs a pseudorandom number generator. That's a deterministic formula that starts from a seed and churns out numbers that pass every statistical test for randomness. For sorting and simulation this is plenty, and it has one lovely bonus: reuse the seed and you can replay a run exactly. But there is a real caveat for the armor story. If an attacker can learn or guess your seed, they can predict every "coin" you'll flip, and the defense evaporates. When randomness must survive a determined adversary, use the secrets module, which draws from the operating system's unpredictable entropy instead.

So how much does a single coin flip actually buy? Let's stop asserting and start counting →

02Randomness is armor — the randomized pivot

Count the fixed-pivot disaster first, because the arithmetic is small enough to check by hand. Partitioning a block of m elements costs m−1 comparisons, one for every element measured against the pivot. If every pivot is the extreme value, which is exactly what a sorted list guarantees, the block sizes march down n−1, n−2, …, 1. Add them up and the total is (n−1)+(n−2)+…+1 = n(n−1)/2, which is the textbook O(n²). And it is not approximate. I ran a comparison-counting quicksort on a sorted list of n = 2000, and the fixed-pivot version reported 1,999,000 comparisons — precisely n(n−1)/2, to the comparison.

Now make the pivot random. Each pivot is equally likely to be any rank, so the partitions come out balanced on average. The expected number of comparisons works out to about 2n ln n (≈ 1.39·n·log₂n), squarely in the O(n log n) class. I fed it the same sorted input, this time with a random pivot: about 24,600 comparisons, averaged over twenty runs. That's roughly 80× fewer comparisons, and it happened on the very input that murdered the deterministic version.

One phrase in that sentence deserves a hard look: on average. An expected value is just a long-run average, the number your results settle toward over many repeats. Roll a fair six-sided die many times and the average roll settles toward 3.5, even though no single roll can be 3.5. Check it: (1+2+3+4+5+6)/6 = 21/6 = 3.5. The randomized quicksort claim is the same kind of statement. Any single run might draw unlucky pivots and land above that 24,600 average, and that's fine. The promise is about the average over the coin flips, and that average holds for every input, sorted or not.

The deep move: we didn't change the input, we changed who chooses. The chance that a random shuffle happens to land on any one specific bad arrangement is 1/n!. For a mere 20 items that's 1 in 2.4×10¹⁸; for a shuffled deck of 52 cards it's 1 in 8×10⁶⁷ — more arrangements than there are atoms in the galaxy. The worst case didn't vanish. It became something you will not witness before the heat death of the universe.

def quicksort_count(a, pivot="first"):
    a = list(a); comps = 0; stack = [(0, len(a)-1)]
    while stack:
        lo, hi = stack.pop()
        if lo >= hi: continue
        if pivot == "random":                 # ← the one line that armors it
            p = random.randint(lo, hi); a[p], a[hi] = a[hi], a[p]
        pv = a[hi]; i = lo
        for j in range(lo, hi):
            comps += 1                         # count every key comparison
            if a[j] < pv:
                a[i], a[j] = a[j], a[i]; i += 1
        a[i], a[hi] = a[hi], a[i]
        stack.append((lo, i-1)); stack.append((i+1, hi))
    return comps
# sorted input, n=2000:  fixed → 1,999,000   random → ~24,600  (≈80× fewer)

Line by line, here is what the code is doing. We sort a copy with an explicit stack instead of recursion, so the sorted worst case can't overflow Python's call stack. The if pivot == "random" line swaps a randomly chosen element into the pivot slot, and that one swap is the entire defense. The inner for loop is the partition. The comps += 1 line tallies each key comparison, so we measure the cost directly instead of trusting a formula. Everything else is ordinary Lomuto partitioning, the same scheme you stepped in T28.

import random
from statistics import mean

random.seed(2024)
N = 52


def fisher_yates(a):
    """Deal from the fan: each round pulls ONE uniform card into a settled slot."""
    a = list(a)
    for i in range(len(a) - 1, 0, -1):
        j = random.randint(0, i)          # 0..i INCLUSIVE — i+1 choices
        a[i], a[j] = a[j], a[i]
    return a


def riffle(deck):
    """One human riffle: cut near the middle, then drop cards from the two
       packets — the fuller packet is likelier to drop next (Gilbert-Shannon-Reeds)."""
    n = len(deck)
    cut = sum(random.random() < 0.5 for _ in range(n))   # binomial cut
    left, right = deck[:cut], deck[cut:]
    out = []
    while left and right:
        if random.random() < len(left) / (len(left) + len(right)):
            out.append(left.pop(0))
        else:
            out.append(right.pop(0))
    return out + left + right


def rising_sequences(perm):
    """How many increasing runs the deck decomposes into. A brand-new deck
       has 1. A truly random deck of n averages (n+1)/2."""
    pos = [0] * len(perm)
    for where, card in enumerate(perm):
        pos[card] = where
    runs = 1
    for c in range(1, len(perm)):
        if pos[c] < pos[c - 1]:
            runs += 1
    return runs


new = list(range(N))
print("brand-new deck        :", rising_sequences(new), "rising sequence")

TRIALS = 4000
for k in (1, 2, 3, 5, 7, 10):
    seqs = []
    for _ in range(TRIALS):
        d = list(new)
        for _ in range(k):
            d = riffle(d)
        seqs.append(rising_sequences(d))
    print(f"{k:2d} riffle(s)          : mean {mean(seqs):6.2f}   max seen {max(seqs):3d}"
          f"   theory caps it at 2**{k} = {2**k}")

fy = [rising_sequences(fisher_yates(new)) for _ in range(TRIALS)]
sh = []
for _ in range(TRIALS):
    d = list(new)
    random.shuffle(d)
    sh.append(rising_sequences(d))
print(f"your fisher_yates     : mean {mean(fy):6.2f}   max seen {max(fy):3d}")
print(f"random.shuffle        : mean {mean(sh):6.2f}   max seen {max(sh):3d}")
print(f"a uniform deck should : mean {(N + 1) / 2:6.2f}")

# fairness of your own shuffle, measured on a deck small enough to count
from collections import Counter
tally = Counter(tuple(fisher_yates("ABC")) for _ in range(600_000))
print(sorted((("".join(k)), v) for k, v in tally.items()))


# ---------- the actual run, python 3.12.7 (seed 2024) ----------
#
# brand-new deck        : 1 rising sequence
#  1 riffle(s)          : mean   2.00   max seen   2   theory caps it at 2**1 = 2
#  2 riffle(s)          : mean   4.00   max seen   4   theory caps it at 2**2 = 4
#  3 riffle(s)          : mean   7.96   max seen   8   theory caps it at 2**3 = 8
#  5 riffle(s)          : mean  19.60   max seen  26   theory caps it at 2**5 = 32
#  7 riffle(s)          : mean  24.70   max seen  32   theory caps it at 2**7 = 128
# 10 riffle(s)          : mean  26.29   max seen  34   theory caps it at 2**10 = 1024
# your fisher_yates     : mean  26.56   max seen  35
# random.shuffle        : mean  26.52   max seen  34
# a uniform deck should : mean  26.50
# [('ABC', 99844), ('ACB', 100036), ('BAC', 99861), ('BCA', 99776), ('CAB', 100285), ('CBA', 100198)]
#
#
# ---------- reading it ----------
#
# 1. Your fisher_yates measured 26.56 and random.shuffle measured 26.52,
#    against a theoretical 26.50. They are the same algorithm -- read
#    inspect.getsource(random.Random.shuffle) and you will find your own
#    loop, running downward with an inclusive roll.
#
# 2. Three riffles: 7.96 rising sequences where a fair deck has 26.5. The
#    deck is still about 70% ordered. It FELT random because your hands
#    were unpredictable; the deck was not, and those are different claims.
#
# 3. Look at the max-seen column against 2**k. One riffle can never produce
#    more than 2 rising sequences -- a riffle interleaves two packets, so
#    it can at most double the runs. Three riffles cap at 8, and the
#    measurement never exceeded 8. That ceiling is not a statistical
#    tendency, it is arithmetic: k riffles cannot manufacture more than
#    2**k rising sequences, however hard you shuffle. This is the Bayer-
#    Diaconis result, and it is why the usual advice is seven riffles.
#
# 4. The cap stops binding around k = 6 (2**6 = 64 is already past what a
#    52-card deck can show), which is why the means keep creeping upward
#    after that -- 24.70 at seven riffles, 26.29 at ten -- instead of
#    jumping. Randomness arrives gradually; the ceiling lifts all at once.
#
# 5. The ABC census: six orderings, 600,000 trials, every count within
#    0.3% of 100,000. That is what "each of the n! orderings with
#    probability 1/n!" looks like when you actually count it.
#
# 6. random.seed(2024) makes every number above reproducible on your
#    machine. Delete the seed and the third digit moves; the story does
#    not.

Armor is one use of a coin. The other is stranger: using randomness to compute a number you can't otherwise reach →

03Compute by throwing darts — Monte Carlo

Now turn the whole idea inside out, because randomness has a second job. So far we used it to protect an exact algorithm from a hostile input. But we can also use it to estimate an answer that's too hard to compute directly. The technique is called Monte Carlo, after the casino, and its motto is worth memorizing: to learn about a system, sample random outcomes from it and count.

The cleanest demonstration computes π. Draw a 1×1 square, and inside it the quarter of a circle of radius 1 centered at one corner. The square's area is 1; the quarter-circle's area is (1/4)·π·r² = π/4. Throw darts uniformly at random into the square. The fraction that land inside the quarter circle must, on average, equal the ratio of the areas — π/4. So π ≈ 4 × (darts inside) / (total darts). A dart at (x, y) is inside exactly when x² + y² ≤ 1 — one multiply, one compare, no trigonometry, and the symbol π appears nowhere in the code. You measure it by counting.

Let's make the counting concrete before we run the code. Suppose you throw 1,000 darts and 787 of them land inside the quarter circle. Your estimate is 4 × 787 / 1000 = 3.148, which is about 0.006 above the true 3.14159…. Notice what did the work there: pure counting, one hit ratio, one multiply by 4. The uniformity of the darts is the load-bearing assumption, because darts that favored one corner would skew the hit ratio. That's why random.random() matters here: it spreads its draws evenly across [0,1), so every patch of the square gets hit in proportion to its area.

def estimate_pi(n):
    inside = 0
    for _ in range(n):
        x, y = random.random(), random.random()   # a dart in the unit square
        if x*x + y*y <= 1.0:                        # inside the quarter circle?
            inside += 1
    return 4 * inside / n                           # fraction inside ≈ π/4

#      n         estimate      error
#     100        3.120000     0.021593
#    1000        3.128000     0.013593
#   10000        3.126000     0.015593   ← note: WORSE than n=1000!
#  100000        3.137280     0.004313
# 1000000        3.140592     0.001001
# 10000000       3.142310     0.000717

Four lines do the whole job. random.random() draws a coordinate in [0,1), so the pair (x, y) is a uniform dart. The one comparison x*x + y*y <= 1.0 asks "is it inside?", inside tallies the hits, and 4 * inside / n turns the ratio into π. The table is a real run, and look closely at it: at n = 10,000 the estimate was actually worse than at n = 1,000. That's not a bug. A single random run can get lucky or unlucky, and no single run owes you accuracy. What improves reliably is never one run — it's the typical error across many, and pinning that down is the next section.

04The √N law — and pricing the future

Here's the law that governs every Monte Carlo estimate, and it's worth burning into memory: the typical error shrinks like 1/√N, where N is the number of samples. I measured it directly. Run the π estimator 200 times at N = 100 and the estimates scatter with a standard deviation of about 0.165. Quadruple to N = 400 and the scatter halves to 0.081. Quadruple again to 1600 → 0.039; again to 6400 → 0.020. Every ×4 in samples halves the error — the fingerprint of 1/√N.

You can see the halving inside the formula itself. Quadruple the samples and the error goes from 1/√N to 1/√(4N), and √(4N) = 2√N. So the new error is exactly half the old one, straight from the square root. Check it against the measurements: half of 0.165 is 0.0825, and the measured scatter at N = 400 was 0.081. Half of 0.081 is 0.0405, and the next measurement came in at 0.039. The law isn't a slogan, because it's sitting right there in the run data.

Flip that around and it stings. To win one more correct decimal digit, meaning a 10× smaller error, you need 100× more samples. Monte Carlo hands you the first digits cheaply and then charges brutally for the last ones. So why 1/√N? Each sample is an independent draw, and averaging N independent draws is exactly the setup where the spread of the average narrows as 1/√N. It's the same statistics that make a poll of 1,000 people accurate to a few percent, not a few thousandths. None of this is special to π — it's the signature of averaging independent randomness, anywhere it appears.

And this is exactly why Wall Street runs on Monte Carlo. The fair price of a financial option is an expected value — the average payoff across every possible future path of the underlying stock, discounted back to today. For the simplest options there's a closed-form answer, the famous Black–Scholes formula. But for the messy real ones, path-dependent, many assets, exotic payoffs, the integral has no formula at all. So you don't solve it — you sample it: simulate thousands of random futures, compute the payoff in each, and average. I priced a European call both ways to see this land, and Black–Scholes says 10.4506. Monte Carlo, averaging simulated futures, gave 11.02 at 1,000 paths, then 10.49 at 100,000 paths, then 10.4458 at a million paths. That's the √N law crawling toward the true price. Risk desks simulate millions of paths every night to answer one question: "how much could we lose tomorrow?" That is expected-value thinking made computational: when you cannot reason about the average, sample it.

If you've never met an option, here's the whole contract in one breath. A European call gives you the right, on one fixed expiry date, to buy a stock at a pre-agreed strike price K. Say the strike is K = 100 and the stock finishes at 112. You exercise, buy at 100, and the option was worth 12. If the stock finishes at 93 instead, you simply walk away, and the option pays 0. That's why the payoff formula is max(ST − K, 0): the final price minus the strike, floored at zero. The hard part is that today nobody knows where the stock will finish, so the fair price has to be an average over all the possible futures.

def mc_call(n, S0, K, r, sigma, T):
    drift = (r - 0.5*sigma**2)*T
    vol   = sigma*(T**0.5)
    total = 0.0
    for _ in range(n):
        z  = random.gauss(0, 1)               # one random shock
        ST = S0*math.exp(drift + vol*z)        # one simulated future price
        total += max(ST - K, 0.0)              # the call's payoff on that path
    return math.exp(-r*T) * total/n            # discounted average payoff
# Black–Scholes = 10.4506;  mc_call(1_000_000) → 10.4458

The loop is the whole idea, so read it slowly. random.gauss(0, 1) draws one random "shock" from a bell curve. Then S0*exp(drift + vol·z) turns that shock into one possible price at expiry. And max(ST − K, 0) is what the call pays on that particular future. Average over a million futures, discount back to today, and you've evaluated an integral no one wrote down — by sampling it.

Quicksort's coins and π's darts gamble in opposite directions. Naming that difference is the whole taxonomy of randomized algorithms →

05Two ways to gamble — Las Vegas vs Monte Carlo

Look at what the coins actually put at risk in each algorithm. Randomized quicksort always returns a correctly sorted list, because the coins only affect how long it takes — a bad streak makes it slower, never wrong. The π estimator is the mirror image: it always finishes in exactly N steps, but the answer is only approximate. So a bad streak there makes it wrong-ish, never slow. Those are the two families of randomized algorithm, and knowing which one you're holding tells you exactly what can go wrong.

The families have names, and fittingly, both names are casinos. The always-correct, sometimes-slow family is called Las Vegas, and randomized quicksort lives there. You always walk out with the right answer, and only the waiting time varies. The fixed-time, approximately-right family is called Monte Carlo, and the π estimator lives there. The runtime is locked, but the answer carries an error you must respect. So when you meet a new randomized algorithm, ask the one sorting question first: did the coins put the time at risk, or the answer?

There's a beautiful bridge between the two families. If you can check an answer cheaply, you can convert Monte Carlo into Las Vegas. Guess with a fast randomized method, then verify the guess, and if it's wrong, simply guess again. The verify step turns "probably right" into "definitely right, eventually" — you've traded a small error probability for a small amount of extra runtime.

06A fair sample from an endless stream — reservoir sampling

The problem sounds impossible when you first hear it. You're watching a firehose: every tweet, every log line, every sensor reading. You want to keep a uniformly random sample of, say, one item. The requirement: when the stream finally ends, every item that ever flowed past had an equal chance of being the one you kept. Now the catch, and it's a triple one. You don't know how long the stream is, because it could be billions of items or never end at all. You get exactly one pass over the data, and you have O(1) memory — room for the one item you're holding, nothing more.

The solution is one line of genius, and it's called Reservoir sampling: keep the first item. When the i-th item arrives, replace whatever you're holding with it, with probability exactly 1/i. That's the entire algorithm, and notice something strange: n, the stream's length, never appears anywhere in the rule.

def reservoir_one(stream):
    kept = None
    for i, item in enumerate(stream, start=1):   # i = 1, 2, 3, ...
        if random.random() < 1.0 / i:            # newest barges in w.p. 1/i
            kept = item
    return kept
# length-10 stream, 1,000,000 trials — how often each item was the survivor:
#   item 0: 10.00%   item 1: 10.04%   item 2:  9.98%   item 3:  9.99%
#   item 4:  9.97%   item 5: 10.00%   item 6: 10.06%   item 7:  9.98%  ...
# every item ≈ 1/10, dead fair — with room for exactly one item in memory

The loop reads one item at a time, holding a single survivor. At position i, the test random.random() < 1.0/i comes up true with probability 1/i. So item 1 is kept for sure, item 2 has a 1-in-2 shot at bumping it, item 3 a 1-in-3 shot, and so on down the stream. The comment in the code is a real million-trial run. Every one of the ten items ended up the survivor almost exactly a tenth of the time, deviating by at most 0.06 percentage points. Perfectly fair, in one pass, in constant memory.

Before the general proof, run the smallest interesting case by hand: a stream of just three items. Item 3 is the survivor only if it's kept at its own step, which happens with probability 1/3. Item 2 must be kept at step 2 and then survive item 3, so its chance is 1/2 × 2/3 = 1/3. Item 1 is kept for sure, then must survive both later steps: 1 × 1/2 × 2/3 = 1/3. All three land on exactly 1/3, with nothing tuned and no stream length known in advance. The general proof below is this same multiplication, written once for every position.

Here's the magic, worked out. Consider item i. It gets kept at step i with probability 1/i. To still be the survivor at the end, it must then not be replaced by items i+1, i+2, …, n. Item j replaces with probability 1/j, so it fails to replace with probability (j−1)/j. Multiply the whole chain of survival:

(1/i) · (i/(i+1)) · ((i+1)/(i+2)) · … · ((n−1)/n) = 1/n

and it telescopes: every numerator cancels the denominator just before it, leaving exactly 1/n. Every item ends with probability 1/n, whether it arrived first or billionth. The math is fair by construction, and the length n is never needed to make it so.

you type
$ python random_idiom.py

import random, secrets, inspect

random.seed(7)                                  # replayable coins

deck = ["A", "K", "Q", "J", "10"]

# ---------- 1. shuffle: in place, returns None ----------
d = list(deck)
print(random.shuffle(d), d)

# ---------- 2. sample: k WITHOUT replacement, a new list ----------
print(random.sample(deck, 3), deck)

# ---------- 3. choices: k WITH replacement, weights allowed ----------
print(random.choices(deck, k=5))
print(random.choices(["hit", "miss"], weights=[1, 9], k=10))

# ---------- 4. choice / randrange / randint / uniform / gauss ----------
print(random.choice(deck), random.randrange(5), random.randint(0, 4))
print(round(random.uniform(0, 1), 4), round(random.gauss(0, 1), 4))

# ---------- 5. the seed makes a run replayable ----------
random.seed(42); a = random.sample(range(100), 5)
random.seed(42); b = random.sample(range(100), 5)
print(a, b, a == b)

# ---------- 6. CPython's shuffle IS Fisher-Yates. read it, don't trust me ----------
print(inspect.getsource(random.Random.shuffle).strip())

# ---------- 7. the two traps ----------
try:
    random.sample({"A", "K", "Q"}, 2)
except TypeError as err:
    print("TypeError:", err)
try:
    random.sample(deck, 9)
except ValueError as err:
    print("ValueError:", err)

# ---------- 8. when an adversary is watching, random is not enough ----------
print(secrets.choice(deck) in deck, len(secrets.token_hex(8)))
sysrand = random.SystemRandom()
print(len(sysrand.sample(range(100), 3)))

you see
None ['10', 'A', 'J', 'K', 'Q']
['A', '10', 'K'] ['A', 'K', 'Q', 'J', '10']
['Q', '10', 'K', 'A', 'Q']
['miss', 'miss', 'hit', 'miss', 'miss', 'miss', 'miss', 'hit', 'miss', 'hit']
K 0 4
0.8585 -0.1375
[81, 14, 3, 94, 35] [81, 14, 3, 94, 35] True
def shuffle(self, x):
        """Shuffle list x in place, and return None."""

        randbelow = self._randbelow
        for i in reversed(range(1, len(x))):
            # pick an element in x[:i+1] with which to exchange x[i]
            j = randbelow(i + 1)
            x[i], x[j] = x[j], x[i]
TypeError: Population must be a sequence.  For dicts or sets, use sorted(d).
ValueError: Sample larger than population or is negative
True 16
3

Step back and see the shape of this whole chapter. Randomness bought us three different things, and the first was armor, because no adversary can aim at a random pivot. The second was computation: sample the darts, the futures, the light rays, and average your way to an answer no formula gives. The third was fairness under ignorance, a uniform sample drawn without ever knowing how much data there is. The transferable move, the one the 1% run on reflex, is this. When being predictable is your weakness, add a coin. When the exact answer is out of reach, estimate it by trying many random cases and averaging. And when you want a property like fairness, work backwards from the property to the rule that produces it. Stanislaw Ulam invented Monte Carlo not by out-computing the neutron equations but by asking a humbler question. Recovering from illness and playing solitaire, he wondered about his odds of winning at the cards. He realized that instead of grinding through the combinatorics, he could just deal many random hands and count. That question — "what if I just try it thousands of times and count?" — is available to you on every hard problem you'll ever meet.

Six chapters of patterns; one volume of thinking. Next we assemble them into a single operating system for the algorithmic mind — how to look at any new problem and know, in seconds, which technique it wants. Chapter 47 is the difference between the 99% and the 1%, and it closes the volume →