Problem 213: Flea Circus

Detailed mathematical approach

Problem Summary

A \(30\times30\) grid of squares holds \(900\) fleas, one on each square. Each time a bell rings, every flea jumps to a square that shares a side with its own, chosen at random: as usual, each of the \(4\) neighbours of an interior square, the \(3\) of an edge square or the \(2\) of a corner is equally likely, and every choice is independent of the earlier choices and of the other fleas. We must find the expected number of unoccupied squares after \(50\) rings, rounded to six decimal places; the statement gives no intermediate values. ^[1]

The whole solution in one sentence: the fleas jump independently, so a square is empty with a probability equal to the product, over all \(900\) fleas, of the probability that this flea is elsewhere, and by linearity of expectation the answer is the sum of these products over all squares, where the probabilities of each flea come from spreading its position over the grid one ring at a time.

Mathematical Approach

Following all fleas together is hopeless, since their joint positions form \(900^{900}\) arrangements, more than \(10^{2658}\). An expectation, however, can be built flea by flea and square by square, with every step checked on small grids.

1. One flea is a random walk

Number the rows and columns from \(0\) to \(29\) and write \(\deg(u)\) for the number of neighbours of square \(u\): it is \(2\) for the \(4\) corners, \(3\) for the \(4\cdot28=112\) other edge squares and \(4\) for the \(28^2=784\) interior squares. A flea on \(u\) moves to each neighbour with probability \(1/\deg(u)\), whatever its earlier path, so it performs a Markov chain, the simple random walk on the grid. ^{[2] [3]} Let \(p_s^{(t)}(v)\) be the probability that the flea starting on square \(s\) is on square \(v\) after \(t\) rings; \(p_s^{(0)}\) is \(1\) on \(s\) and \(0\) elsewhere. To be on \(v\) after \(t+1\) rings, the flea must be on a neighbour \(u\sim v\) after \(t\) rings and then choose \(v\):

$$p_s^{(t+1)}(v)=\sum_{u\sim v}\frac{p_s^{(t)}(u)}{\deg(u)}.$$

Every square splits its probability into equal shares, one for each neighbour. Nothing is lost, so the values of one flea always add up to \(1\), and the values after \(t\) rings are the entries of the \(t\)-th power of the transition matrix. ^[2] On a \(3\times3\) grid, a flea starting in a corner is after one ring on each of the two edge squares beside it with probability \(1/2\). These pass \(1/6\) to each of their three neighbours, so after two rings the flea is back in its corner with probability \(1/6+1/6=1/3\), in the centre with probability \(1/3\) and in each of the two nearest corners with probability \(1/6\).

2. One square: a product over the fleas

Let \(q(v)\) be the probability that square \(v\) is empty after \(50\) rings, that is, that every flea is elsewhere; in the small examples below, \(q(v)\) refers to their own number of rings. Whether the flea from \(s\) is elsewhere depends only on the choices of that flea, so these \(900\) events are mutually independent, and the probability that all of them happen is the product of their probabilities: ^[4]

$$q(v)=\prod_{s}\bigl(1-p_s^{(50)}(v)\bigr).$$

The product needs each \(p_s^{(50)}(v)\) on its own, not only their sum \(\lambda(v)\), the expected number of fleas on \(v\). Two fleas that each reach \(v\) with probability \(1/2\) leave it empty with probability \(1/4\), while a flea that surely comes and one that never comes leave it empty with probability \(0\), although \(\lambda(v)=1\) in both cases. That is why the programs follow the \(900\) starting squares one at a time.

3. All squares: linearity of expectation

Let \(I_v\) be \(1\) if square \(v\) is empty and \(0\) otherwise, so that the number of empty squares is \(X=\sum_v I_v\). The expectation of an indicator is the probability of its event, and the expectation of a finite sum of random variables with finite expectations, such as these indicators, is the sum of their expectations, whether they are independent or not. ^[5] Therefore

$$\begin{aligned}\mathbb E[X]&=\sum_{v}\mathbb E[I_v]\\&=\sum_{v}\prod_{s}\bigl(1-p_s^{(50)}(v)\bigr),\end{aligned}$$

which is what the programs evaluate. The formula never asks how the emptiness of one square affects another, and that effect can be strong.

4. Worked example: the \(2\times2\) grid

Call the squares \(\mathtt A\), \(\mathtt B\) in the top row and \(\mathtt C\), \(\mathtt D\) in the bottom row; each square has the two squares of the other diagonal as its neighbours. After one ring the fleas from \(\mathtt B\) and \(\mathtt C\) are each on \(\mathtt A\) with probability \(1/2\), and the fleas from \(\mathtt A\) and \(\mathtt D\) cannot be there, so

$$q(\mathtt A)=\Bigl(1-\frac12\Bigr)^2(1-0)^2=\frac14.$$

By symmetry \(\mathbb E[X]=4\cdot\tfrac14=1\), the checkpoint of the C++ program. Yet \(\mathtt A\) is empty only if the fleas from \(\mathtt B\) and \(\mathtt C\) both jump to \(\mathtt D\), so \(\mathtt A\) and \(\mathtt D\) are never empty together, although \(\tfrac14\cdot\tfrac14\ne0\); linearity does not mind. Of the \(2^4=16\) equally likely combinations of jumps, \(4\), \(8\) and \(4\) leave \(0\), \(1\) and \(2\) squares empty, with mean \(1\). After any number \(t\ge1\) of rings each flea is on the two squares of one diagonal with probability \(1/2\) each, so every square is again reached by two fleas with probability \(1/2\) each, and \(\mathbb E[X]=1\) for every \(t\ge1\).

5. Worked example: the \(3\times3\) grid after one ring

After one ring a flea can be on \(v\) only if it started on a neighbour \(u\) of \(v\), and then it is there with probability \(1/\deg(u)\). So \(q(v)\) is the product of \(1-1/\deg(u)\) over the neighbours \(u\) of \(v\):

• a corner has two edge neighbours with three neighbours each, so \(q=(2/3)^2=4/9\);
• the middle square of an edge has two corners and the centre as neighbours, so \(q=\tfrac12\cdot\tfrac12\cdot\tfrac34=\tfrac3{16}\);
• the centre has four edge neighbours, so \(q=(2/3)^4=16/81\).

With four squares of each of the first two kinds,

$$\begin{aligned}\mathbb E[X]&=4\cdot\frac49+4\cdot\frac3{16}+\frac{16}{81}\\&=\frac{883}{324}\approx2.725309.\end{aligned}$$

Listing all \(2^4\cdot3^4\cdot4=5184\) equally likely combinations of jumps of the nine fleas gives exactly \(883/324\) as well. On the real grid the same one-ring rule gives \(162769/576\approx282.585069\).

6. The chessboard colouring

Colour the squares like a chessboard, with \((r,c)\) black when \(r+c\) is even. Every jump changes \(r+c\) by \(1\) and so the colour: after an even number of rings a flea stands on its starting colour, and the walk is said to have period two. ^[3] Hence \(p_s^{(50)}(v)=0\) whenever \(s\) and \(v\) have different colours, so only the \(450\) fleas of the colour of \(v\) can reach it, and the other \(450\) factors of \(q(v)\) equal \(1\). Distance matters too, since a flea needs at least \(|r-r'|+|c-c'|\) jumps from \((r,c)\) to \((r',c')\), and opposite corners are \(58\) jumps apart. Of the \(405\,000\) pairs \((s,v)\) of the same colour, \(520\) are more than \(50\) jumps apart; every other pair is joined by a walk of exactly \(50\) jumps, a shortest path followed by jumps back and forth. So exactly \(404\,480\) of the \(810\,000\) values \(p_s^{(50)}(v)\) are positive.

7. Corners and centre after 50 rings

Summing the rule of Section 1 over all starting squares shows that the expected numbers of fleas \(\lambda^{(t)}(v)=\sum_s p_s^{(t)}(v)\) obey the same rule, started from \(1\) on every square. Without a border, every square would always receive a quarter of a flea from each of its four neighbours, and \(\lambda\) would stay \(1\). The border breaks this balance. If a flea's position has the distribution \(\pi(v)=\deg(v)/3480\), where \(3480\) is the sum of all degrees, each neighbour \(u\) of \(v\) sends \(\pi(u)/\deg(u)=1/3480\) to \(v\), so \(v\) again receives \(\pi(v)\); this distribution is therefore stationary, and it gives a corner only half the weight of an interior square. ^[3] Stationarity says nothing about \(50\) rings, but the rule of Section 1 does: after \(50\) rings a corner holds \(0.542556\) fleas on average and is empty with probability \(0.579342\), while \((14,14)\), one of the four central squares, holds \(1.002120\) fleas and is empty with probability \(0.365932\).

8. A shortcut that is not accurate enough

Since \(1-p\le e^{-p}\) for every \(p\), each square satisfies \(q(v)\le e^{-\lambda(v)}\), the probability that a Poisson variable with mean \(\lambda(v)\) is zero. As \(\ln(1-p)=-p-p^2/2-\dots\), the product is smaller by a factor of about \(e^{-S/2}\) with \(S=\sum_s p_s^{(50)}(v)^2\). For the square \((14,14)\), \(e^{-\lambda}=0.367101\), \(S=0.006340\) and \(e^{-\lambda-S/2}=0.365939\), but the exact product is \(0.365932\): the first approximation is wrong in the third decimal and even the corrected one in the sixth, so the products must be evaluated as they stand.

9. Exact fractions and floating point

Every jump divides a probability by \(2\), \(3\) or \(4\), all divisors of \(12\), so the scaled values \(a_s^{(t)}(v)=12^t\,p_s^{(t)}(v)\) start as \(0\) or \(1\) and follow the integer rule

$$a_s^{(t+1)}(v)=\sum_{u\sim v}\frac{12}{\deg(u)}\,a_s^{(t)}(u).$$

With \(D=12^{50}\), each \(q(v)\) is the product of the \(900\) integers \(D-a_s^{(50)}(v)\) divided by \(D^{900}\), so the expected value is an exact fraction. The programs use floating point instead: in binary64 with rounding to nearest, an operation whose nonzero result lies in the normal range, as all nonzero results here do, has a relative error of at most \(2^{-53}\approx1.1\times10^{-16}\). ^[6] All values are nonnegative and the only subtraction, \(1-p\), has \(p\lt0.044\), so no cancellation magnifies these errors. The exact fraction, computed with Python integers, shows that the double-precision loops of the programs, summed from left to right, are off by \(1.2\times10^{-13}\), mostly in the final sum. The exact value is more than \(4.8\times10^{-7}\) away from any point where its rounding to six decimals would change, so the printed decimals are those of the exact value.

How the Code Works

The C++ program ^[7] reads --size= and --steps=, with the defaults \(30\) and \(50\), through parse_int_after_prefix and accepts --skip-checkpoints; an unknown argument or a size of \(0\) gives exit code \(1\). A size of \(1\) is accepted although a lone square has no neighbour; from one ring on, the division by zero makes the flea's probability vanish and the result is \(1\), and very large values overflow the unchecked int arithmetic. Unless skipped, run_checkpoints requires solve(2, 1) to be within \(10^{-15}\) of \(1\), the value of Section 4, or the exit code is \(2\). Then main prints the result as a double with six decimals.

solve numbers square \((r,c)\) as r * size + c and keeps the neighbours of each square, above, below, left and right, in neighbors, a vector of vectors whose lengths are the degrees. The vector empty_prob starts at \(1\) everywhere and collects the products \(q(v)\). For each start, dist is \(1\) on start and \(0\) elsewhere. Each ring zeroes next, adds dist[i] divided by the number of neighbours of i to every neighbour of each square i, and swaps the two vectors without copying; pushing shares out of each square forms the same sums as collecting them in Section 1. After the last ring each entry of empty_prob is multiplied by 1 - dist[i], and after all \(900\) starts the sum of the entries is \(\mathbb E[X]\).

Squares with p == 0.0L, by Section 6 at least half of them at every ring, are skipped. The test is exact: an unreachable square never receives a share, and a reachable one has at least the probability \(4^{-50}\approx7.9\times10^{-31}\) of a single path, far above the smallest positive double.

C++ stores the probabilities as long double, which is at least as precise as double: the 80-bit extended format with GCC on x86 Linux, but plain double with Microsoft Visual C++ and on 64-bit ARM under macOS, ^[8] and by Section 9 that is enough. Python's solve and Java's main run the same loops in the same order with binary64 numbers, for the fixed values \(30\) and \(50\) and without checkpoint. Python creates a new list nxt at every ring, and since version 3.12 its sum adds floats more accurately, ^[9] which can change the last bits of the total. Java swaps its double[] arrays through tmp, and its printf follows the default locale, so a German or Turkish setting gives a decimal comma. ^[10]

Complexity and Verification

For each of the \(N=n^2\) starting squares the program runs \(T\) rings over the \(N\) squares with at most four additions each, so for \(T\ge1\) the time is \(O(TN^2)=O(Tn^4)\), and the memory is \(O(N)\). Here the inner loop runs \(40\,500\,000\) times; \(12\,362\,684\) of these passes find a nonzero probability and make \(48\,045\,096\) additions, and \(810\,000\) multiplications form the products. On our machine the programs take about \(0.03\) seconds in C++, \(0.1\) in Java and \(2\) in Python. Repeated squaring of the transition matrix would need seven products of \(900\times900\) matrices with \(900^3\approx7.3\times10^8\) multiplications each, nearly a hundred times more work.

Besides the C++ checkpoint, we ran the following checks, which are not part of the published programs.

• Listing all equally likely combinations of jumps gives \(1\) on the \(2\times2\) grid after \(1\), \(2\) and \(3\) rings and \(883/324\) on the \(3\times3\) grid after one ring.
• The exact joint distribution of the set of occupied squares, built flea by flea, uses independence but neither the product nor linearity; it matches the formula on the \(3\times3\) grid after \(1\), \(2\) and \(3\) rings and on the \(4\times4\) grid after \(1\) and \(2\) rings, for example \(44759/15552\approx2.878022\) for \(3\times3\) and two rings, which --size=3 --steps=2 also prints.
• \(10^6\) simulated runs on a \(10\times10\) grid with \(10\) rings give \(36.1361\pm0.0031\) empty squares, against \(36.135630\) from the formula, and for the real problem \(10^5\) runs differ from the formula by \(1.25\) standard errors. With a run-to-run standard deviation of about \(9.35\), bringing the standard error \(9.35/\sqrt M\) of \(M\) runs down to \(5\times10^{-7}\) would take about \(3.5\times10^{14}\) runs.
• The three programs print the six decimals of the exact fraction of Section 9 (C++ with Clang and GCC, Java 17, Python 3.11, 3.13 and 3.14); with Python 3.11, and with C++ on our ARM machine, where long double is double, the totals agree bit for bit.

Footnotes and References

The references below supply the background facts; the worked examples, the counts and the checks were computed for this article.

1. Project Euler 213 — Flea Circus. The official statement: 900 fleas on a 30 × 30 grid, one per square, each jumping at every ring of a bell to an adjacent square at random; it asks for the expected number of unoccupied squares after 50 rings, rounded to six decimal places.
2. Wikipedia — Markov chain. The next state depends only on the present state, and the k-th power of the transition matrix, whose rows sum to one, gives the k-step transition probabilities.
3. D. A. Levin and Y. Peres — Markov Chains and Mixing Times, second edition. Section 1.4 defines the simple random walk on a graph, a step to a uniformly chosen neighbour; Section 1.3 treats chains of period two, which alternate between two classes of states; Example 1.12 shows that deg(y)/2|E| is a stationary distribution of the walk.
4. Wikipedia — Independence (probability theory). For mutually independent events, the probability that all of them happen is the product of their probabilities.
5. Wikipedia — Expected value. Expectation is linear, E[X + Y] = E[X] + E[Y] whenever the right side is well defined, without any independence assumption; the expectation of an indicator is the probability of its event.
6. Wikipedia — Machine epsilon. For binary64 with rounding to nearest, the unit roundoff, or rounding machine epsilon, 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶ bounds the relative rounding error; the gap between 1 and the next larger number is twice as large.
7. C++ — Euler213.cpp. The linked immutable C++ revision contains parse_arguments, solve and run_checkpoints. It is the implementation record for this article.
8. Wikipedia — Long double. It only has to be at least as precise as double; it is 80-bit extended with GCC on x86 Linux, but the same as double with Microsoft Visual C++ and on 64-bit ARM under Windows and macOS.
9. Python tutorial — Floating-Point Arithmetic. Almost all platforms map Python floats to IEEE 754 binary64 values with 53 bits of precision; the documentation of the built-in sum notes that floats are summed more accurately since Python 3.12.
10. The Java Language Specification, Java SE 17 — §4.2. The values of double are exactly the 64-bit IEEE 754 binary64 numbers; PrintStream.format, which printf calls, always uses the default locale of the FORMAT category.