Let's flip a coin

The last chapter ended with some odd advice. To find the weakest spot in a network, instead of computing flows you can merge random lines until only two cities are left, and repeat that many times. It sounds like giving up on thinking. Yet a coin can be the best tool a programmer has. It gives answers where exact calculation is hopeless, it protects against inputs crafted by an attacker, and it lets you keep track of millions of things in a few kilobytes. This chapter is laid out like a casino. Each table has its own game, and at every one the house wins. The house is us.

Table one: solitaire

That is how the Monte Carlo method got its name: to find a probability, play many times and count the share of successes. We met it in Chapter 12, where Kristen Nygaard played out the fates of neutrons for the first Norwegian reactor. And Ulam turned up in Chapter 22 with a question about increasing subsequences of a random permutation; he was always asking how random things behave.

Canfield is a bit much to start with: it has many rules and choices of moves, so the probability depends on the player as well. We’ll take Clock patience instead, a game with no choices at all. The deck is dealt face down into 13 piles of 4 cards: twelve piles in a circle like a clock face, for the aces, twos and so on up to the queens, and a thirteenth in the middle for the kings. Turn up the top card of the middle pile and slip it under the pile of its rank: a seven goes to seven o’clock, say. From there take the top face-down card, put it under its own pile, and so on. The game comes out if all 52 cards get turned up. Working out the exact probability of that is a job for combinatorics. We’ll do what Ulam did.

The game ends when we come to an empty pile, and only the middle pile can be the empty one. We arrive at any other pile only with a card of its rank, and there are four such cards, as many as the face-down cards in the pile. The middle pile, though, needs one card more: we started there without arriving with a king. So the game always breaks off at the fourth king, and the only question is whether any cards are still face down at that moment.

Here the exact answer is known, $1/13$, and it can be proved. But the machine took about a second to arrive at roughly 0.077, without a single formula. For Canfield, where a lot also depends on the player’s decisions, no simple formula is in sight at all, and dealing is all that is left. Below is the same experiment on our server, with a picture. The widget’s second game, Frustration, is simpler still: you turn the cards up one at a time while calling out “ace, two, three, …, king, ace, …” and lose as soon as a card matches the rank you call. Peter Doyle, Charles Grinstead and Laurie Snell worked out its odds of winning: about one game in 61.6.

How many deals it takes

Each deal is a toss of a biased coin that comes up “success” with probability $p$. After $n$ tosses the share of successes averages $p$, and its spread, the standard deviation, is $\sqrt{p(1-p)/n}$. By the central limit theorem, with probability about 95% the share lands within two such deviations of $p$. For Clock with a thousand deals that is $2\sqrt{\frac{1}{13} \cdot \frac{12}{13} / 1000} \approx 0.017$, so an answer of “somewhere from 0.06 to 0.094” is on the rough side. With a hundred thousand deals it is $\pm 0.0017$.

That is slow. But the $1/\sqrt{n}$ law has a strong side too: it contains neither the difficulty of the problem nor its number of dimensions. Solitaire with 52 cards or with a thousand, a neutron in three dimensions or the price of an option that depends on a hundred stocks: the error shrinks the same way, like $1/\sqrt{n}$, and only the factor in front, the spread of a single trial, depends on the problem. We’ll check this on a number everybody knows.

Table two: darts

Hang up a square target with side 1 and draw on it a quarter of a circle of radius 1 centered at one corner. The quarter circle has area $\pi/4$, the square has area 1. If darts are thrown at random, so that every point of the square is equally likely, the share of throws that land in the quarter circle will be close to $\pi/4$. Multiply the share by four, and you have $\pi$. A point $(x, y)$ lies in the quarter circle if $x^2 + y^2 \le 1$.

A million darts, and only two or three correct digits after the decimal point. For $\pi$ that is laughable: formulas give trillions of digits. Our math course, “Mathematics, the Queen of the Sciences,” has a funnier story, about Buffon’s needle and the Italian Lazzarini, who in 1901 got six correct digits “by random throws,” having chosen the moment to stop suspiciously well.

Nobody needs Monte Carlo for $\pi$; the method earns its keep on problems that have no formula. You could also find the area of a shape with a grid: test the centers of $\sqrt n$ by $\sqrt n$ cells, and for a smooth shape in the plane the error will be about $1/\sqrt n$, no worse. But in ten dimensions a grid of a million points has fewer than four points along each axis: $4^{10} \approx 10^6$. Such a grid sees next to nothing. Random points still give an error of $1/\sqrt n$. That is why Monte Carlo computes volumes in spaces of many dimensions, lighting in computer graphics (every ray of light is a random history, like a neutron) and the risks of financial portfolios. In the task “The area of a blot” at the end of the chapter, you’ll compute the area of a shape for which a formula is awkward.

Table three: roulette

Solitaire and darts give an approximate answer. But a coin also helps algorithms that must answer exactly. Here there are two ways of taking a risk, and each has its own name.

A Las Vegas algorithm always answers correctly, and only its running time depends on the coin. A Monte Carlo algorithm runs in a predictable time but may be wrong, with a small probability that we control. The name “Las Vegas” was coined in 1979 by the Hungarian mathematician László Babai as a partner for the already familiar “Monte Carlo”: one casino against another.

The random quicksort of Chapter 21 is a Las Vegas algorithm. It always sorts the list, and the random choice of pivot affects only the number of comparisons. There we saw that with the first element as the pivot, a sorted list costs quadratic time, and we promised to explain where random quicksort’s average of about $1.39\,n\log_2 n$ comparisons comes from. Time to settle that debt.

The proof says nothing about the input: the average is taken over the coin tosses, not over lists. An attacker who knows our code can’t craft a bad input, because they don’t know which pivots we will pick. It is the same protection as the salt in the hash table of Chapter 16. We’ll check the theorem by counting. The exact average is known: $2(n+1)H_n - 4n$, where $H_n = 1 + \frac12 + \ldots + \frac1n$.

All five runs land within a few percent of the expected value. The ratio to $n \log_2 n$ is still below 1.39, about 1.2 at a hundred thousand, because $2n\ln n$ has a correction of roughly $2.8n$ taken off it, and at modest $n$ that correction shows. The spread is small, and one can prove that the probability of making even twice the average number of comparisons quickly goes to zero as $n$ grows.

We have already met an example of a Monte Carlo algorithm, at the end of the last chapter. Contracting random edges finds a minimum cut of a graph with $n$ vertices with probability at least $\frac{2}{n(n-1)}$. That isn’t much, but repeat the experiment $n^2$ times and keep the best cut, and the chance of failure falls below $e^{-2}$; with $n^2 \ln n$ repetitions, below $1/n^2$. Repetition helps any Monte Carlo algorithm: if each run fails with probability $q$, independently of the others, then all $k$ runs fail together with probability $q^k$.

Table four: dice

This table tests whether a number is prime. For small numbers, trying divisors up to $\sqrt n$ is enough, as in the task from Chapter 5. But an RSA key, the subject of Chapter 60, is built from primes hundreds of digits long, and a three-hundred-digit number has $10^{150}$ candidates below its square root. Trial division is out.

The first idea came from Fermat’s little theorem, which you’ll find in the math course: if $n$ is prime, then $a^{n-1} \equiv 1 \pmod n$ for every $a$ not divisible by $n$. Modular exponentiation is fast: the function pow(a, e, n) does it in $O(\log e)$ multiplications. So we pick a random $a$ and look: if $a^{n-1} \not\equiv 1$, then $n$ is certainly composite. Such an $a$ is called a witness to compositeness. The trouble is that some composite numbers have almost no Fermat witnesses: the Carmichael numbers. The smallest of them is $561 = 3 \cdot 11 \cdot 17$.

With $a = 2$, Fermat’s test sees a 1 and is ready to call 561 prime. But the chain holds evidence. We wrote $n - 1 = 560$ as $2^4 \cdot 35$ and raised 2 first to the odd power 35, then squared four times; the last number of the chain is $2^{560}$. Right before the first 1 stands 67: $67^2 \equiv 1 \pmod{561}$, although $67 \not\equiv \pm 1$. Modulo a prime this can’t happen. If $x^2 \equiv 1 \pmod p$, then $p$ divides $(x-1)(x+1)$, and a prime that divides a product divides one of the factors, so $x \equiv 1$ or $x \equiv -1$. Having found an “extra” square root of 1, we have caught 561 red-handed. The evidence even gives away the divisors: $\text{gcd}(67 - 1, 561) = 33$ and $\text{gcd}(67 + 1, 561) = 17$.

This is what the Miller–Rabin test is built on. Write $n - 1 = 2^s d$ with $d$ odd and compute the chain $a^d, a^{2d}, \ldots, a^{2^s d}$. If $n$ is prime, the chain either starts with 1 straight away or meets $-1$ (that is, $n - 1$) somewhere, and after that come only ones. Anything else is evidence. In 1976 Gary Miller proposed a deterministic test whose correctness was proved on the assumption of the unproved extended Riemann hypothesis. In 1980 Michael Rabin turned it into a probabilistic test that needs no hypotheses at all. Rabin, and independently Louis Monier, proved what the test rests on: for an odd composite number, at most a quarter of the numbers $a$ from 1 to $n - 1$ “lie,” passing the check without being witnesses. The other three quarters are witnesses.

Miller–Rabin is a Monte Carlo algorithm with one-sided error. The answer “composite” is always right, since a witness is a proof. The answer “prime” may be false, but with twenty random values of $a$ the probability of that is at most $4^{-20} < 10^{-12}$. The number $2^{521} - 1$ was checked in milliseconds. In January 1952 Raphael Robinson proved it prime on the SWAC computer with a different test, the Lucas–Lehmer test, which works only for numbers of the form $2^p - 1$. It was the first Mersenne prime found by a computer.

Primes for keys are found this way: take a random odd number of the right length and keep testing until a prime turns up. Primes are common enough (near $N$, roughly one number in $\ln N$ is prime) that for three-hundred-digit numbers a few hundred tries are enough on average. In 2002 Manindra Agrawal, Neeraj Kayal and Nitin Saxena found a deterministic test that runs in polynomial time with no hypotheses at all. It was a breakthrough in theory, but in practice people still roll the dice: the probabilistic test is much faster.

Table five: the lottery drum

Tickets ride past the table on a conveyor belt, and a winner has to be drawn so that every ticket has the same chance. Nobody knows how many tickets there will be, the belt can’t be rewound, and your hand holds one ticket. Everyday tasks look like this too: pick a random line from a hundred-gigabyte server log, a post from a live feed, a word from a novel without loading the novel into memory.

The solution fits in three lines. Hold one ticket in your hand. When the $i$-th ticket comes by, swap yours for it with probability $1/i$. The first ticket we take for certain, the second with probability ½, the third with probability ⅓, and so on.

The six letters in the program’s last line were each chosen about ten thousand times, as the theorem promises. To choose $k$ items, hold $k$ tickets: take the first $k$ outright, and let the $i$-th ticket replace a random one in your hand with probability $k/i$. This is reservoir sampling. In Knuth’s book it is called “Algorithm R,” and Knuth credits it to Alan Waterman; in 1985 Jeffrey Vitter analyzed the whole family of such algorithms and found a way to skip many items at once instead of tossing a coin for each one. In the task “The reservoir” you’ll build one with $k$ places.

Table six: the bouncer

At the casino door stands a bouncer with a list of ten million unwanted guests. He can’t remember every face, so he has been allowed to make mistakes, but only in one direction: now and then he may turn away an ordinary guest, but he must never let in an unwanted one. A browser that checks an address against a list of dangerous sites has the same problem, and so does a database that wants to know, before an expensive read from disk, whether the key is there at all. The set from Chapter 8 never errs, but it stores the items themselves, tens of bytes apiece.

In 1970 Burton Bloom described in Communications of the ACM a structure that gets by on about ten bits per item. His example was hyphenation. Of half a million English words, nine in ten can be hyphenated by simple rules, and the rest need a dictionary on a slow disk. If memory holds a compact “list of exceptions” that sometimes errs on the side of “yes, this is an exception,” then a needless trip to the disk is rare, and nothing is ever missed.

A Bloom filter is an array of $m$ bits, all zero at first, plus $k$ different hash functions, each of which maps an item to a bit position. To add an item, set all its $k$ bits to one. To test an item, inspect the same $k$ bits: if even one is zero, the item was certainly never added; if all are ones, it “probably was.” Those ones may have been set by other items, and then the filter is wrong: this is a false positive. The opposite mistake can’t happen, because bits are only ever set, never cleared.

How many hash functions should we take? Too few, and each word is checked by one or two bits, so a stranger slips through easily. Too many, and the filter quickly fills up with ones. Here is the count. After $n$ additions, a given bit is still zero if none of the $kn$ hash values landed on it, which has probability $(1 - 1/m)^{kn} \approx e^{-kn/m}$. A stranger gets through if all its $k$ bits are ones:

The minimum comes at $k = \frac{m}{n}\ln 2$, when ones fill half the bits, and then $p \approx 0.6185^{\,m/n}$. For one percent of errors you need about 9.6 bits per item, for one in a thousand about 14.4, and this depends neither on the number of items nor on their length. We’ll test it on War and Peace: the filter remembers the words of the first half of the novel, and we ask it about words that occur only in the second half.

Experiment and formula agree, and the minimum falls at $k = 7$, as promised: $10 \ln 2 \approx 6.9$. The method positions doesn’t compute seven different hashes: it takes two numbers from a single SHA-256 hash and builds the rest from them as $h_1 + i\,h_2$. Adam Kirsch and Michael Mitzenmacher proved in 2006 that this shortcut hardly hurts the filter. The special method __contains__, one of those with two underscores from Chapter 12, lets us write w in bloom, as with an ordinary set. And we spend a whole byte per bit only to keep things simple: a working filter packs eight bits into every byte, and Chapter 28 will show how that is done.

Bloom filters work inside the databases Bigtable, Cassandra and HBase: before reading a file from disk, the database asks the filter whether the key could be in it. The content delivery network Akamai kept a filter of the addresses that had already been requested at least once and cached a file only on its second request, so pages that are viewed only once never made it into the cache. You can’t delete from an ordinary filter, because a bit may have been set by another item as well; filters with counters in place of bits were invented for that.

Table seven: the head count

The last table is by the door. How many different people came in during the day, if many of them came in several times? For a website this is the number of unique visitors, for a search engine the number of distinct queries, for a linguist the number of distinct words in a book. The exact answer means remembering everyone you have seen, a set as large as the answer. Can we get by with a kilobyte?

Hashing again. Run every visitor through a good hash function: it turns each of them into a random-looking string of bits, and the same visitor always gives the same string, so repeat visits change nothing. Remember a single number: the largest count of zeros that any string has started with. A hash that starts with ten zeros turns up in about one of every $2^{10} = 1024$ different visitors. If we have seen one, there were probably on the order of a thousand different visitors.

A single such number is a very rough estimate: one zero more or less doubles it. So the visitors are spread over $m$ cells by the first bits of their hash, each cell remembers its own maximum, and then the maxima are carefully averaged. Philippe Flajolet and Nigel Martin proposed a similar “probabilistic counting” in 1985, and in 2007 Flajolet and his coauthors Éric Fusy, Olivier Gandouet and Frédéric Meunier refined it into the HyperLogLog algorithm, with an error of about $1.04/\sqrt{m}$. Each cell holds a number up to 64, which takes 6 bits. In the Redis database a HyperLogLog keeps 16,384 cells: 12 kilobytes and an error of about 0.8%. Many analytics systems count unique visitors the same way.

A thousand cells, less than a kilobyte, and the number of different words in War and Peace comes out with an error of under two percent. A set of all the words would take megabytes. The operations >> and & cut the bits we need out of a number; Chapter 28 covers them in detail. The constant alpha and the correction for empty cells come from the 2007 paper. Averaging the maxima overstates the answer a little, and the constant was chosen to make up for that.

Cashing out

Time to count the winnings. At every table the coin gave us something that exact calculation would not, and for every win we paid a clear price.

Now a confession: there was never any coin. The random module is a pseudorandom number generator, a deterministic program that turns a starting number, the “seed,” into a long sequence that looks random. In Python it is the Mersenne Twister, published by Makoto Matsumoto and Takuji Nishimura in 1998. The same seed gives the same sequence, which is handy for experiments: the result can be reproduced.

Tasks

Four tasks from four of the seven tables. In two of them the answer is random, and the tests check it statistically: many runs, and the shares have to land where the theory says they should.

What next

Michael Rabin, whose name is on the primality test, turned randomness on one more problem seven years later, together with Richard Karp: finding a pattern in a text. The problem sounds modest until you see the scale. The human genome is about three billion letters from the alphabet A, C, G, T, and a biologist needs to find a short motif in it, and then thousands of others. The brute-force way, laying the pattern against every position in turn, makes billions of comparisons for each pattern. An automaton that never steps back, a rolling hash and a tree of letters will do it faster. That is the subject of the next chapter.