Problem 213: Flea Circus

Mathematical Development

Definition 1. Let $G = (V, E)$ be the grid graph with $V = \{(r, c) : 0 \leq r, c < 30\}$ and edges between horizontally or vertically adjacent squares. Let $N = |V| = 900$. Index the squares $1, \ldots, N$.

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 vertex $i$ and $j \sim i$ means $j$ is adjacent to $i$.

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

Proof. For any vertex $i$, $\sum_{j} T_{ij} = \sum_{j \sim i} 1/\deg(i) = \deg(i) \cdot (1/\deg(i)) = 1$. $\square$

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

Proof. By induction on $k$. The base case $k = 1$ holds by definition of $T$. For the inductive step, $(T^{k+1})_{ij} = \sum_{\ell} (T^k)_{i\ell} \cdot T_{\ell j}$, which by the inductive hypothesis equals $\sum_{\ell} P(\text{at } \ell \text{ after } k \text{ steps from } i) \cdot P(\text{move to } j \text{ from } \ell) = P(\text{at } j \text{ after } k+1 \text{ steps from } i)$, by the law of total probability. $\square$

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

Proof. Step 1: By linearity of expectation (which holds without any independence assumption),

Step 2: Square $i$ is empty if and only if none of the $N$ fleas occupies it. Since fleas move independently (each flea’s random walk is independent of all others), the events $\{\text{flea } j \neq i \text{ after 50 rounds}\}$ are mutually independent across distinct $j$. Therefore:

where $(T^{50})_{ji}$ is the probability that flea $j$ (starting at square $j$) is at square $i$ after 50 steps, by Theorem 1. $\square$

Remark. To avoid numerical underflow in the product $\prod_j (1 - (T^{50})_{ji})$, one may compute $\sum_j \log(1 - (T^{50})_{ji})$ and exponentiate. The function $\log(1 - x)$ is well-conditioned for small $x$, and the library function log1p(-x) provides enhanced precision.

Editorial

We build transition matrix T. We then compute T^R via matrix exponentiation. Finally, accumulate expected empty squares. We enumerate the admissible parameter range, discard candidates that violate the derived bounds or arithmetic constraints, and update the final set or total whenever a candidate passes the acceptance test.

Pseudocode

Build transition matrix T
Compute T^R via matrix exponentiation
Accumulate expected empty squares

Complexity Analysis

Time. Matrix exponentiation via repeated squaring requires $O(\log R)$ matrix multiplications, each costing $O(N^3)$. Thus the total is $O(N^3 \log R)$. For $N = 900$ and $R = 50$, this is approximately $900^3 \cdot 6 \approx 4.4 \times 10^9$. Alternatively, propagating each of the $N$ flea distributions individually over $R$ rounds costs $O(N \cdot R \cdot N) = O(N^2 R) = O(900^2 \cdot 50) \approx 4 \times 10^7$, which is more practical. The product computation adds $O(N^2)$. Overall: $O(N^2 R)$.

Space. $O(N^2)$ for storing the full probability matrix, or $O(N)$ if accumulating products column-by-column.

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()