Problem 213: Flea Circus

Mathematical Development

Definition 1. Let $G = (V, E)$ be the grid graph with

and edges between horizontally or vertically adjacent squares. Let $N = |V| = 900$.

Definition 2. The transition matrix $T \in \mathbb{R}^{N \times N}$ of the random walk on $G$ is defined by

where $\deg(i) \in \{2,3,4\}$ is the degree of square $i$.

Lemma 1 (Row stochasticity). The matrix $T$ is row-stochastic: $\sum_j T_{ij} = 1$ for all $i$.

Proof. For a fixed square $i$, exactly $\deg(i)$ entries in row $i$ are nonzero, each equal to $1/\deg(i)$. Their sum is therefore 1. $\square$

Theorem 1 (Chapman-Kolmogorov). The entry $(T^k)_{ij}$ equals the probability that a random walk starting at square $i$ is at square $j$ after exactly $k$ steps.

Proof. The case $k = 1$ is the definition of $T$. The inductive step is the usual matrix-product form of the law of total probability. $\square$

Theorem 2 (Expected empty squares). Let $\mathbf{1}_i$ be the indicator that square $i$ is empty after 50 rounds. Then

Proof. By linearity of expectation,

Square $i$ is empty exactly when none of the 900 fleas lands there. Since the fleas move independently,

using Theorem 1 for the individual landing probabilities. $\square$

Remark. When multiplying many factors of the form $1-p$, a logarithmic accumulation

can be numerically more stable than repeated floating-point multiplication.

Editorial

The empty-square expectation separates cleanly by square. For one fixed target square, each flea independently misses it with some probability, so the probability that the square ends up empty is a product over all starting squares. Summing those 900 emptiness probabilities gives the final expectation.

The real work is therefore to compute the 50-step distribution of a flea starting from each square. That can be done either by matrix powers or, more practically, by propagating one probability vector for 50 rounds from each starting square. After each source distribution is computed, multiply every board position by the probability that this particular flea is absent there. At the end, those products are exactly the emptiness probabilities.

Pseudocode

Label the 900 squares and precompute the legal neighbors of each one.
Initialize absent_probability[i] = 1 for every square i.

For each starting square s:
    Set a probability vector current with current[s] = 1.

    Repeat 50 times:
        Create a zero vector next.
        For each square u:
            Split current[u] equally among the neighbors of u
            and add those contributions into next.
        Replace current by next.

    For each target square i:
        absent_probability[i] *= 1 - current[i]

answer = sum of absent_probability[i] over all squares i
Return answer

Complexity Analysis

• Time: Propagating one flea distribution for 50 rounds costs $O(50N)$ state updates, so handling all 900 starting squares costs $O(N^2 R)$ with $N = 900$ and $R = 50$.
• Space: $O(N)$ for the current and next probability vectors, plus $O(N)$ for the accumulated absence probabilities.

Answer

#include <bits/stdc++.h>
using namespace std;

int main() {
    const int R = 30, C = 30;
    const int N = R * C;
    const int ROUNDS = 50;

    auto idx = [&](int r, int c) { return r * C + c; };

    int dr[] = {-1, 1, 0, 0};
    int dc[] = {0, 0, -1, 1};

    vector<int> deg(N);
    for (int r = 0; r < R; r++)
        for (int c = 0; c < C; c++) {
            int d = 0;
            for (int k = 0; k < 4; k++)
                if (r + dr[k] >= 0 && r + dr[k] < R && c + dc[k] >= 0 && c + dc[k] < C)
                    d++;
            deg[idx(r, c)] = d;
        }

    vector<double> prod_not(N, 1.0);

    for (int sr = 0; sr < R; sr++) {
        for (int sc = 0; sc < C; sc++) {
            vector<double> cur(N, 0.0), nxt(N, 0.0);
            cur[idx(sr, sc)] = 1.0;

            for (int t = 0; t < ROUNDS; t++) {
                fill(nxt.begin(), nxt.end(), 0.0);
                for (int r = 0; r < R; r++)
                    for (int c = 0; c < C; c++) {
                        int i = idx(r, c);
                        if (cur[i] == 0.0) continue;
                        double p = cur[i] / deg[i];
                        for (int k = 0; k < 4; k++) {
                            int nr = r + dr[k], nc = c + dc[k];
                            if (nr >= 0 && nr < R && nc >= 0 && nc < C)
                                nxt[idx(nr, nc)] += p;
                        }
                    }
                swap(cur, nxt);
            }

            for (int i = 0; i < N; i++)
                prod_not[i] *= (1.0 - cur[i]);
        }
    }

    double answer = 0.0;
    for (int i = 0; i < N; i++)
        answer += prod_not[i];

    cout << fixed << setprecision(6) << answer << endl;
    return 0;
}

import numpy as np

def solve():
    R, C = 30, 30
    N = R * C
    ROUNDS = 50

    def idx(r, c):
        return r * C + c

    dr = [-1, 1, 0, 0]
    dc = [0, 0, -1, 1]

    deg = np.zeros(N, dtype=float)
    for r in range(R):
        for c in range(C):
            deg[idx(r, c)] = sum(
                1 for k in range(4) if 0 <= r + dr[k] < R and 0 <= c + dc[k] < C
            )

    T = np.zeros((N, N), dtype=float)
    for r in range(R):
        for c in range(C):
            i = idx(r, c)
            for k in range(4):
                nr, nc = r + dr[k], c + dc[k]
                if 0 <= nr < R and 0 <= nc < C:
                    T[i, idx(nr, nc)] = 1.0 / deg[i]

    Tn = np.linalg.matrix_power(T, ROUNDS)
    log_prod_not = np.sum(np.log1p(-Tn), axis=0)
    answer = np.sum(np.exp(log_prod_not))

    print(f"{answer:.6f}")

if __name__ == "__main__":
    solve()