Problem 280: Ant and Seeds

Project Euler Problem 280 Solution

EulerSolve provides an optimized solution for Project Euler Problem 280, Ant and Seeds, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary An ant performs a uniform random walk on a $5\times5$ grid. Initially the five bottom-row cells each contain one seed, the five top-row cells are empty, and the ant starts at the center $(2,2)$. If the ant is not carrying a seed and reaches a bottom cell that still contains one, it picks it up immediately. If it is carrying a seed and reaches an empty top cell, it drops the seed immediately. We want the expected number of steps until all five seeds have been transported to the top row. Mathematical Approach 1) What information actually matters? A naive Markov state would record the ant position together with the exact locations of all five seeds. That is far more detailed than necessary, because seeds only matter when they sit in the bottom row waiting to be picked up, or in the top row after being delivered. So a state is completely described by: $$B\subseteq\{0,1,2,3,4\}\quad\text{: columns whose bottom seed is still present},$$ $$T\subseteq\{0,1,2,3,4\}\quad\text{: columns already occupied in the top row},$$ the ant position $p\in\{0,\dots,24\}$, and a carrying flag. The code stores $B$ and $T$ as 5-bit masks. This is already a major reduction: there are only $25$ possible positions and $32$ masks per row....

Detailed mathematical approach

Problem Summary

An ant walks on a $5\times5$ grid, starting on the central square. At every step it moves to a square that shares a side with its present square, chosen with equal probability among the $2$, $3$ or $4$ such squares. At the start each square of the bottom row holds a seed. When the ant carries no seed and reaches a bottom square that still holds one, it picks the seed up; it drops it on the first empty square of the top row that it reaches. We must find the expected number of steps until all five seeds have been dropped in the top row, rounded to $6$ decimal places; no smaller cases are given. ^[1]

The whole solution in one sentence: between two pick-ups or drops the ant makes an ordinary random walk until it first reaches one of a few target squares, so the mean length of each such trip and the probability of each end square come from small linear systems, and a table over the seeds still at the bottom and the filled top squares chains the ten trips together.

Mathematical Approach

We fix what the process must remember, compute the length and the end of a single trip, checked by hand on a $3\times3$ grid, and then join the trips.

1. Squares and states

As in the code, the columns $x=0,\dots,4$ run from left to right and the rows $y=0,\dots,4$ from top to bottom; square $(x,y)$ has the number $p=5y+x$, so the ant starts on square $12$. We write $t_j=j$ and $b_j=20+j$ for the top and the bottom square of column $j$, and $d(u)$ for the number of neighbours of square $u$: $2$ for the $4$ corners, $3$ for the other $12$ border squares and $4$ for the $9$ inner squares.

The future of the walk depends on four things: the ant's square, whether it carries a seed, the set $B$ of columns whose bottom square still holds its seed, and the set $T$ of columns whose top square holds a seed. Every seed is at the bottom, carried, or at the top for good, and the ant carries at most one, so $|B|+|T|=5$ without a seed and $|B|+|T|=4$ with one. By Vandermonde's identity ^[2] there are

$$\sum_{k=0}^{5}\binom5k\binom5{5-k}=\binom{10}5=252$$

pairs $(B,T)$ of the first kind and $\binom{10}4=210$ of the second, so with the $25$ squares there are $6300+5250=11550$ combinations. The $1260$ of them with the ant on a seed it must pick up, or carrying on an empty top square, never occur, since events happen on arrival; on the other $10290$ the walk is a Markov chain that stops once $T=\{0,1,2,3,4\}$. Its expected number of steps solves a linear system with one unknown per state; ^[3] the verification below does that, but trips make the work far smaller.

2. Trips between events

Call a pick-up or a drop an event; it takes no step. Without a seed the ant ignores top squares and empty bottom squares, so the next event happens the first time it stands on some $b_j$ with $j\in B$. With a seed it cannot take another one and ignores full top squares, so the next event happens the first time it stands on some $t_j$ with $j\notin T$. A trip runs from the start or from an event to the next event, so it is an ordinary random walk until it first enters a fixed nonempty target set $S$. Five pick-ups and five drops make ten trips, the first one from the centre; for each we need its expected length and the probability of each end square.

3. The length of a trip: a $3\times3$ example

Take a $3\times3$ grid with its bottom row as the target, and let $h(u)$ be the expected number of steps from $u$ until the walk first stands on the bottom row, so $h=0$ there. From any other square the walk steps to a random neighbour $v$ and then needs $h(v)$ more steps on average, so $h(u)$ is $1$ plus the average of $h$ over the neighbours of $u$. By left–right symmetry there are four unknowns: $\alpha$ at the top corners, $\beta$ at the top middle square, $\gamma$ at the middle squares of the two sides and $\delta$ at the centre. The equations are

$$\begin{aligned}\alpha&=1+\tfrac12(\beta+\gamma),\\ \beta&=1+\tfrac13(2\alpha+\delta),\\ \gamma&=1+\tfrac13(\alpha+\delta+0),\\ \delta&=1+\tfrac14(2\gamma+\beta+0),\end{aligned}$$

where the zeros come from neighbours on the bottom row. Subtracting the third equation from the second gives $\beta=\gamma+\alpha/3$; then the first equation becomes $\alpha=\tfrac65(1+\gamma)$, the fourth becomes $\delta=\tfrac1{20}(22+17\gamma)$, and the third turns into $19\gamma=106$. Hence

$$\begin{aligned}\alpha&=\tfrac{150}{19},&\beta&=\tfrac{156}{19},\\ \gamma&=\tfrac{106}{19},&\delta&=\tfrac{111}{19},\end{aligned}$$

which satisfy all four equations. From the centre the ant needs on average $111/19\approx5.84$ steps to reach the bottom row.

4. The general rule, and why it has one solution

For any nonempty target set $S$, this conditioning on the first step, called first-step analysis, gives ^[4]

$$h_S(u)=\begin{cases}0,&u\in S,\\1+\dfrac1{d(u)}\displaystyle\sum_{v\sim u}h_S(v),&u\notin S,\end{cases}$$

where $v\sim u$ means that $v$ is a neighbour of $u$. On the $5\times5$ grid every trip ends, with finite expected length: every square is at most $4+4=8$ steps from $S$, and the walk follows a given path of at most $8$ steps with probability at least $4^{-8}$, so a trip lasts more than $8k$ steps with probability at most $(1-4^{-8})^k$, and its mean is at most $8\cdot4^8$. The system also has only one solution. The difference $g$ of two solutions is $0$ on $S$ and, outside $S$, equals the average of its values at the neighbours. If its largest value $M$ is taken outside $S$, all neighbours of that square have the value $M$ too, since an average equals its largest term only if all terms are equal; following a path into $S$, we find $M$ on $S$, where $g=0$. So $M=0$, likewise the smallest value is $0$, and $g=0$. This is the maximum principle for such averaging, or harmonic, functions; ^[5] equivalently, $I-Q$ is invertible, $Q$ being the moves between squares outside $S$. ^[3]

5. Where a trip ends

For $s\in S$, let $q_{S,s}(u)$ be the probability that a trip from $u$ first enters $S$ at $s$. The first step gives ^[4]

$$q_{S,s}(u)=\begin{cases}1,&u=s,\\0,&u\in S,\ u\ne s,\\\dfrac1{d(u)}\displaystyle\sum_{v\sim u}q_{S,s}(v),&u\notin S.\end{cases}$$

The matrix is the one of Section 4 and only the right-hand side differs, so one elimination yields $h_S$ and every $q_{S,s}$; since every trip ends, the $q_{S,s}(u)$ add up to $1$ over $s\in S$. On the $3\times3$ grid with $s$ the middle bottom square, unknowns at the same four places satisfy $\alpha'=\tfrac12(\beta'+\gamma')$, $\beta'=\tfrac13(2\alpha'+\delta')$, $\gamma'=\tfrac13(\alpha'+\delta'+0)$ and $\delta'=\tfrac14(2\gamma'+\beta'+1)$, with the solution $\alpha'=6/19$, $\beta'=7/19$, $\gamma'=5/19$, $\delta'=9/19$. So from the centre the first bottom square reached is the middle one with probability $9/19$ and each corner with probability $\tfrac12(1-9/19)=5/19$.

On the real grid the first trip, from the centre to the bottom row, takes $1137582/48301\approx23.55$ steps on average and ends in the columns $0$ to $4$ with probabilities of about $0.169$, $0.202$, $0.257$, $0.202$, $0.169$.

6. Joining the trips

Let $\tau$ be the length of a trip from square $p$. Whether $\tau=n$ is decided by the first $n$ steps, so $\tau$ is a stopping time, and by the strong Markov property the walk after $\tau$, given that it stands on $s$, is a new random walk from $s$, independent of the trip; ^[6] the seeds after the event are determined by $s$ too. Splitting the steps into $\tau$ and the rest and averaging over the end square gives: the expected number of remaining steps is the expected trip length plus the sum over $s\in S$ of $q_{S,s}(p)$ times the expected number of remaining steps after an event at $s$.

7. The two recurrences

For a set $M$ of columns, let $h^{\downarrow}_M$ and $q^{\downarrow}_{M,j}$ be the expected length and the end probabilities of a trip whose targets are the bottom squares $b_j$, $j\in M$, and let $h^{\uparrow}_M$ and $q^{\uparrow}_{M,j}$ be the same for the top squares $t_j$, $j\in M$. Let $E_0(B,T,p)$ and $E_1(B,T,p)$ be the expected numbers of remaining steps when the ant stands on $p$ without and with a seed, and let $\overline T$ be the set of columns whose top square is empty. Section 6 gives

$$\begin{aligned}E_0(B,T,p)&=h^{\downarrow}_{B}(p)\\&\quad+\sum_{j\in B}q^{\downarrow}_{B,j}(p)\,E_1(B\setminus\{j\},T,b_j),\end{aligned}$$

$$\begin{aligned}E_1(B,T,p)&=h^{\uparrow}_{\overline T}(p)\\&\quad+\sum_{j\in\overline T}q^{\uparrow}_{\overline T,j}(p)\,E_0(B,T\cup\{j\},t_j).\end{aligned}$$

The walk ends with the last drop, so $E_0(\varnothing,\{0,\dots,4\},p)=0$ for every $p$, and the answer is $E_0(\{0,\dots,4\},\varnothing,12)$. With one seed left each sum has a single term with probability $1$: for the last seed in column $s$ and the last empty top square in column $u$, that is $T=\{0,\dots,4\}\setminus\{u\}$,

$$\begin{aligned}E_0(\{s\},T,p)&=h^{\downarrow}_{\{s\}}(p)+h^{\uparrow}_{\{u\}}(b_s),\\E_1(\varnothing,T,p)&=h^{\uparrow}_{\{u\}}(p).\end{aligned}$$

8. An order without cycles

With $c=1$ while the ant carries a seed and $c=0$ otherwise, $r=2|B|+c$ counts the events still to come, and each event lowers it by $1$. The recurrences therefore refer only to values with one event fewer, and starting from the end states, where $r=0$, every value is computed once, in the order $r=1,2,\dots,10$, without iteration. As $|B|+|T|+c=5$, $r=10-2|T|-c$: running through $|T|=4,3,2,1,0$ and, for each, through the states with a seed before those without is this order. Only $2\cdot31=62$ kinds of trips occur, since the targets depend only on $B$ or only on $\overline T$, both nonempty.

9. Worked example: three seeds on the $3\times3$ grid

Now let $t_j$ and $b_j$ be the top and bottom squares of column $j=0,1,2$ of the $3\times3$ grid, with three seeds and the start at the centre. The first trip is the one of Sections 3 and 5, so with $X=E_1(\{1,2\},\varnothing,b_0)$ and $Y=E_1(\{0,2\},\varnothing,b_1)$, and since a pick-up at $b_2$ is the mirror image of one at $b_0$,

$$E=\frac{111}{19}+\frac{10}{19}X+\frac9{19}Y.$$

After a pick-up at $b_1$ the seed goes from the bottom middle square to the top row, the trip from the top middle square to the bottom row turned upside down: it takes $\beta=156/19$ steps on average and ends at $t_1$ with probability $\beta'=7/19$ and at $t_0$ or $t_2$ with probability $6/19$ each. After a drop at $t_1$ the rest can be followed by hand: small systems of the same kind give $19/2$ steps from $t_1$ to a bottom corner and $9$ from there to a top corner, and the last seed then needs one trip between the two corners of a column, $15$ steps on average, and one between opposite corners, $18$ steps, in some order. So $E_0(\{0,2\},\{1\},t_1)=103/2$, the sum of these means. After a drop at $t_0$ or $t_2$ the table gives $4266/85$, so

$$\begin{aligned}Y&=\frac{156}{19}+\frac{7}{19}\cdot\frac{103}{2}+\frac{12}{19}\cdot\frac{4266}{85}\\&=\frac{190189}{3230}\approx58.882.\end{aligned}$$

In the same way $X=1583844/27455\approx57.689$, so $E=66870807/1043290\approx64.096$.

How the Code Works

The C++ program ^[8] numbers the squares as in Section 1 (make_pos, top_pos, bottom_pos), and build_grid stores the neighbours of every square and $d(u)$. A set of columns is a 5-bit mask with bit $j$ for column $j$; $31$ is the full set.

solve_hitting_data(grid, row, mask, ok) writes the systems of Sections 4 and 5 for the targets of mask on row $0$ or $4$ as one $25\times25$ matrix with six right-hand sides: a target square gets an identity row, $0$ for $h$ and $1$ for its own column; another square $u$ gets $1$ on the diagonal, $-1/d(u)$ per neighbour, $1$ for $h$ and $0$ for the probabilities. solve_linear_system applies Gauss–Jordan elimination with partial pivoting: ^[7] per column it swaps the row with the largest remaining absolute entry into place, divides it by the pivot and subtracts multiples of it from the other rows, on all six right-hand sides at once. It skips rows whose factor is below its threshold and would fail if no pivot reached $10^{-18}$; the smallest pivot is about $0.2$. The results go into HittingData records, expect[p] and prob[p][j].

build_solution_data solves the $62$ systems hit_top[mask] and hit_bottom[mask] (masks $1$ to $31$) and fills expected_nc[B][T][p] and expected_c[B][T][p], that is $E_0$ and $E_1$, for the $252+210$ valid mask pairs. From the end values expected_nc[0][31][p] $=0$, top_count runs from $4$ down to $0$ (Section 8), first over the carrying states, whose bottom masks have 4 - top_count bits, with the drop data hit_top[31 ^ top_mask], then over the non-carrying states, whose bottom masks have 5 - top_count bits, with the pick-up data hit_bottom[bottom_mask]. All $25$ squares are computed, although the recurrences and the final lookup read only $b_j$, $t_j$ and the centre; this is cheap and lets the checks cover every square. main prints expected_nc[31][0][12] with six decimals.

Unless --skip-checkpoints is given, run_checkpoints verifies within $10^{-10}$ that the end probabilities of every system add up to $1$ and that $h=0$ on targets, and within $5\cdot10^{-10}$ the one-seed formulas of Section 7 for all $s$, $u$, $p$; a failure ends the program with an error code. solve_expected_steps repeats these steps but is never called. C++ keeps its tables in std::array objects of long double, on some platforms wider than double, and prints the result as a double. Python (nested lists of floats) and Java (double arrays) mirror the functions with the pivot threshold $10^{-15}$; they keep the failure check of the elimination but omit the checkpoints, and Java formats with Locale.US.

Complexity and Verification

Pivot column $k$ changes $25-k$ matrix entries and $6$ right-hand sides in each of the other $24$ rows, so a system needs at most $24\cdot475=11400$ multiply–subtract operations; skipping zero entries brings the total for all $62$ systems down to $347044$. The table has $5250$ carrying and $6275$ non-carrying values besides the end values, each a sum of at most five products, and the two arrays hold $32\cdot32\cdot25=25600$ numbers each. On an Apple M4 Max (clang++ -O2, Python 3.11, OpenJDK 17) the C++ program runs in about $3$ ms, Python and Java in under $0.05$ s including start-up.

All three programs print the same value, and the C++ checks pass. Independent checks for this article, not part of the published solutions:

Sections 4 to 8 in exact rational arithmetic give the same six decimals. The exact value is a fraction with a $135$-digit denominator; the double result of the Python program differs from it by about $5\cdot10^{-12}$, and it lies more than $2\cdot10^{-7}$ from the nearest rounding boundary.
The step-by-step chain of Section 1, whose $10265$ states before the end all occur, solved with a sparse solver, gives the same six decimals.
$10^7$ simulated walks agree with the exact value within $0.4$ standard errors (standard error $0.04$).
For the $3\times3$ version, the trip method and the step-by-step chain with its $246$ states both give exactly $66870807/1043290$, and $10^7$ simulated walks give $64.0956$ with standard error $0.0076$.

Footnotes and References

The references below supply the background facts used above. The trip method, its order of evaluation and all numbers are derived or computed in this article.

Project Euler 280 — Ant and Seeds. The official statement: a random walk on a 5 × 5 grid from the centre; seeds on the bottom row are picked up when the ant carries none and dropped on the first empty top square; the task asks for the expected number of steps until all are in the top row, rounded to 6 decimal places.
Wikipedia — Vandermonde's identity. C(m + n, r) is the sum over k of C(m, k) · C(n, r − k), proved by counting committees of r people taken from m men and n women by their number of men.
Wikipedia — Absorbing Markov chain. With Q the transitions between transient states and R those into absorbing states, N = I + Q + Q² + ⋯ = (I − Q)⁻¹ gives the expected steps before absorption, t = N1, and the absorption probabilities, B = NR.
J. R. Norris, Markov Chains, §1.3 — Hitting times and absorption probabilities. Theorems 1.3.5 and 1.3.2 derive, by conditioning on the first step, the linear equations for the mean hitting times and the hitting probabilities of a set A; Section 5 uses the walk stopped on S, with A = {s}.
P. G. Doyle and J. L. Snell — Random Walks and Electric Networks (arXiv:math/0001057). Sections 1.2.1–1.2.2 show that the probability of reaching one part of the boundary before another is harmonic, an average over the neighbours, and prove the Maximum and Uniqueness Principles.
J. R. Norris, Markov Chains, §1.4 — Strong Markov property. A first hitting time is a stopping time (Example 1.4.1), and by Theorem 1.4.2, given T < ∞ and X_T = i, the chain after a stopping time T starts afresh from i, independently of X₀, …, X_T.
Wikipedia — Gaussian elimination. Gauss–Jordan elimination brings a matrix to reduced row echelon form, and choosing a pivot of largest absolute value improves numerical stability.
C++ — Euler280.cpp. The linked immutable revision is the implementation record for this article; the file shown on this page counts set bits with __builtin_popcount instead of std::popcount.

Problem 280 source code

C++

#include <array>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <string>

namespace {

using i64 = std::int64_t;
using ld = long double;

constexpr int kSize = 5;
constexpr int kCells = kSize * kSize;
constexpr int kCols = 5;
constexpr int kMaskCount = 1 << kCols;
constexpr int kFullMask = kMaskCount - 1;
constexpr int kRhsCount = 1 + kCols;  // expected time + hit probabilities per column

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

int make_pos(const int x, const int y) {
    return y * kSize + x;
}

int top_pos(const int col) {
    return make_pos(col, 0);
}

int bottom_pos(const int col) {
    return make_pos(col, kSize - 1);
}

struct Grid {
    std::array<std::array<int, 4>, kCells> neighbors{};
    std::array<int, kCells> degree{};
};

Grid build_grid() {
    Grid grid;
    for (int y = 0; y < kSize; ++y) {
        for (int x = 0; x < kSize; ++x) {
            const int p = make_pos(x, y);
            int d = 0;
            if (x > 0) {
                grid.neighbors[p][d++] = make_pos(x - 1, y);
            }
            if (x + 1 < kSize) {
                grid.neighbors[p][d++] = make_pos(x + 1, y);
            }
            if (y > 0) {
                grid.neighbors[p][d++] = make_pos(x, y - 1);
            }
            if (y + 1 < kSize) {
                grid.neighbors[p][d++] = make_pos(x, y + 1);
            }
            grid.degree[p] = d;
        }
    }
    return grid;
}

bool solve_linear_system(
    std::array<std::array<ld, kCells>, kCells>& a,
    std::array<std::array<ld, kRhsCount>, kCells>& rhs
) {
    constexpr ld kEps = 1e-18L;

    for (int col = 0; col < kCells; ++col) {
        int pivot = col;
        ld best = std::fabsl(a[col][col]);
        for (int row = col + 1; row < kCells; ++row) {
            const ld cand = std::fabsl(a[row][col]);
            if (cand > best) {
                best = cand;
                pivot = row;
            }
        }
        if (best < kEps) {
            return false;
        }

        if (pivot != col) {
            std::swap(a[pivot], a[col]);
            std::swap(rhs[pivot], rhs[col]);
        }

        const ld inv_pivot = 1.0L / a[col][col];
        for (int j = col; j < kCells; ++j) {
            a[col][j] *= inv_pivot;
        }
        for (int r = 0; r < kRhsCount; ++r) {
            rhs[col][r] *= inv_pivot;
        }

        for (int row = 0; row < kCells; ++row) {
            if (row == col) {
                continue;
            }
            const ld factor = a[row][col];
            if (std::fabsl(factor) < kEps) {
                continue;
            }
            for (int j = col; j < kCells; ++j) {
                a[row][j] -= factor * a[col][j];
            }
            for (int r = 0; r < kRhsCount; ++r) {
                rhs[row][r] -= factor * rhs[col][r];
            }
        }
    }

    return true;
}

struct HittingData {
    std::array<ld, kCells> expect{};
    std::array<std::array<ld, kCols>, kCells> prob{};
};

HittingData solve_hitting_data(const Grid& grid, const int row, const int mask, bool& ok) {
    HittingData data;

    std::array<bool, kCells> is_target{};
    for (int col = 0; col < kCols; ++col) {
        if ((mask & (1 << col)) != 0) {
            is_target[static_cast<std::size_t>(make_pos(col, row))] = true;
        }
    }

    std::array<std::array<ld, kCells>, kCells> a{};
    std::array<std::array<ld, kRhsCount>, kCells> rhs{};

    for (int state = 0; state < kCells; ++state) {
        if (is_target[static_cast<std::size_t>(state)]) {
            a[static_cast<std::size_t>(state)][static_cast<std::size_t>(state)] = 1.0L;
            rhs[static_cast<std::size_t>(state)][0] = 0.0L;
            const int col = state % kSize;
            rhs[static_cast<std::size_t>(state)][static_cast<std::size_t>(1 + col)] = 1.0L;
            continue;
        }

        a[static_cast<std::size_t>(state)][static_cast<std::size_t>(state)] = 1.0L;
        const int deg = grid.degree[static_cast<std::size_t>(state)];
        const ld step = 1.0L / static_cast<ld>(deg);
        for (int i = 0; i < deg; ++i) {
            const int nb = grid.neighbors[static_cast<std::size_t>(state)][static_cast<std::size_t>(i)];
            a[static_cast<std::size_t>(state)][static_cast<std::size_t>(nb)] -= step;
        }
        rhs[static_cast<std::size_t>(state)][0] = 1.0L;
    }

    if (!solve_linear_system(a, rhs)) {
        ok = false;
        return data;
    }

    for (int state = 0; state < kCells; ++state) {
        data.expect[static_cast<std::size_t>(state)] = rhs[static_cast<std::size_t>(state)][0];
        for (int col = 0; col < kCols; ++col) {
            data.prob[static_cast<std::size_t>(state)][static_cast<std::size_t>(col)] =
                rhs[static_cast<std::size_t>(state)][static_cast<std::size_t>(1 + col)];
        }
    }

    return data;
}

struct SolutionData {
    std::array<HittingData, kMaskCount> hit_top{};
    std::array<HittingData, kMaskCount> hit_bottom{};
    std::array<std::array<std::array<ld, kCells>, kMaskCount>, kMaskCount> expected_nc{};
    std::array<std::array<std::array<ld, kCells>, kMaskCount>, kMaskCount> expected_c{};
};

bool build_solution_data(SolutionData& out) {
    const Grid grid = build_grid();

    bool ok = true;
    for (int mask = 1; mask < kMaskCount; ++mask) {
        out.hit_top[static_cast<std::size_t>(mask)] = solve_hitting_data(grid, 0, mask, ok);
        if (!ok) {
            return false;
        }
        out.hit_bottom[static_cast<std::size_t>(mask)] = solve_hitting_data(grid, kSize - 1, mask, ok);
        if (!ok) {
            return false;
        }
    }

    std::array<int, kMaskCount> popcount{};
    for (int mask = 0; mask < kMaskCount; ++mask) {
        popcount[static_cast<std::size_t>(mask)] = __builtin_popcount(static_cast<unsigned>(mask));
    }

    for (int pos = 0; pos < kCells; ++pos) {
        out.expected_nc[0][kFullMask][static_cast<std::size_t>(pos)] = 0.0L;
    }

    for (int top_count = 4; top_count >= 0; --top_count) {
        // Carrying states at this top_count use non-carrying states at top_count + 1.
        for (int top_mask = 0; top_mask < kMaskCount; ++top_mask) {
            if (popcount[static_cast<std::size_t>(top_mask)] != top_count) {
                continue;
            }
            const int empty_top_mask = kFullMask ^ top_mask;
            const HittingData& drop = out.hit_top[static_cast<std::size_t>(empty_top_mask)];

            for (int bottom_mask = 0; bottom_mask < kMaskCount; ++bottom_mask) {
                if (popcount[static_cast<std::size_t>(bottom_mask)] != 4 - top_count) {
                    continue;
                }

                for (int pos = 0; pos < kCells; ++pos) {
                    ld value = drop.expect[static_cast<std::size_t>(pos)];

                    for (int col = 0; col < kCols; ++col) {
                        const int bit = 1 << col;
                        if ((empty_top_mask & bit) == 0) {
                            continue;
                        }
                        value +=
                            drop.prob[static_cast<std::size_t>(pos)][static_cast<std::size_t>(col)] *
                            out.expected_nc[static_cast<std::size_t>(bottom_mask)]
                                           [static_cast<std::size_t>(top_mask | bit)]
                                           [static_cast<std::size_t>(top_pos(col))];
                    }

                    out.expected_c[static_cast<std::size_t>(bottom_mask)]
                                  [static_cast<std::size_t>(top_mask)]
                                  [static_cast<std::size_t>(pos)] = value;
                }
            }
        }

        // Non-carrying states at this top_count use carrying states at same top_count.
        for (int top_mask = 0; top_mask < kMaskCount; ++top_mask) {
            if (popcount[static_cast<std::size_t>(top_mask)] != top_count) {
                continue;
            }
            for (int bottom_mask = 0; bottom_mask < kMaskCount; ++bottom_mask) {
                if (popcount[static_cast<std::size_t>(bottom_mask)] != 5 - top_count) {
                    continue;
                }

                const HittingData& pick = out.hit_bottom[static_cast<std::size_t>(bottom_mask)];
                for (int pos = 0; pos < kCells; ++pos) {
                    ld value = pick.expect[static_cast<std::size_t>(pos)];

                    for (int col = 0; col < kCols; ++col) {
                        const int bit = 1 << col;
                        if ((bottom_mask & bit) == 0) {
                            continue;
                        }
                        value +=
                            pick.prob[static_cast<std::size_t>(pos)][static_cast<std::size_t>(col)] *
                            out.expected_c[static_cast<std::size_t>(bottom_mask ^ bit)]
                                          [static_cast<std::size_t>(top_mask)]
                                          [static_cast<std::size_t>(bottom_pos(col))];
                    }

                    out.expected_nc[static_cast<std::size_t>(bottom_mask)]
                                   [static_cast<std::size_t>(top_mask)]
                                   [static_cast<std::size_t>(pos)] = value;
                }
            }
        }
    }

    return true;
}

bool run_checkpoints(const SolutionData& data) {
    constexpr ld kTol = 1e-10L;

    // Check hitting-system invariants.
    for (int mask = 1; mask < kMaskCount; ++mask) {
        for (const auto* table : {&data.hit_top, &data.hit_bottom}) {
            const HittingData& hit = (*table)[static_cast<std::size_t>(mask)];
            for (int pos = 0; pos < kCells; ++pos) {
                const int row = (table == &data.hit_top) ? 0 : (kSize - 1);
                const int col = pos % kSize;
                const bool is_target = (pos / kSize == row) && ((mask & (1 << col)) != 0);

                ld prob_sum = 0.0L;
                for (int c = 0; c < kCols; ++c) {
                    if ((mask & (1 << c)) != 0) {
                        prob_sum += hit.prob[static_cast<std::size_t>(pos)][static_cast<std::size_t>(c)];
                    }
                }

                if (std::fabsl(prob_sum - 1.0L) > kTol) {
                    std::cerr << "Probability sum checkpoint failed for mask=" << mask
                              << ", pos=" << pos << '\n';
                    return false;
                }
                if (is_target && std::fabsl(hit.expect[static_cast<std::size_t>(pos)]) > kTol) {
                    std::cerr << "Target expectation checkpoint failed for mask=" << mask
                              << ", pos=" << pos << '\n';
                    return false;
                }
            }
        }
    }

    // With one seed left and one top slot empty, the process decomposes into
    // deterministic pick/drop targets.
    for (int seed_col = 0; seed_col < kCols; ++seed_col) {
        const int bottom_mask = 1 << seed_col;
        const int bottom_start = bottom_pos(seed_col);

        for (int empty_col = 0; empty_col < kCols; ++empty_col) {
            const int empty_mask = 1 << empty_col;
            const int top_mask = kFullMask ^ empty_mask;

            const ld carry_from_bottom =
                data.hit_top[static_cast<std::size_t>(empty_mask)]
                    .expect[static_cast<std::size_t>(bottom_start)];

            for (int pos = 0; pos < kCells; ++pos) {
                const ld expected_nc =
                    data.hit_bottom[static_cast<std::size_t>(bottom_mask)]
                        .expect[static_cast<std::size_t>(pos)] +
                    carry_from_bottom;

                const ld actual_nc =
                    data.expected_nc[static_cast<std::size_t>(bottom_mask)]
                                    [static_cast<std::size_t>(top_mask)]
                                    [static_cast<std::size_t>(pos)];
                if (std::fabsl(actual_nc - expected_nc) > 5e-10L) {
                    std::cerr << "One-seed non-carrying checkpoint failed at pos=" << pos
                              << ", seed_col=" << seed_col << ", empty_col=" << empty_col << '\n';
                    return false;
                }

                const ld expected_carry =
                    data.hit_top[static_cast<std::size_t>(empty_mask)]
                        .expect[static_cast<std::size_t>(pos)];
                const ld actual_carry =
                    data.expected_c[0][static_cast<std::size_t>(top_mask)]
                                   [static_cast<std::size_t>(pos)];
                if (std::fabsl(actual_carry - expected_carry) > 5e-10L) {
                    std::cerr << "One-seed carrying checkpoint failed at pos=" << pos
                              << ", empty_col=" << empty_col << '\n';
                    return false;
                }
            }
        }
    }

    return true;
}

ld solve_expected_steps() {
    SolutionData data;
    if (!build_solution_data(data)) {
        return -1.0L;
    }

    if (!run_checkpoints(data)) {
        return -2.0L;
    }

    const int initial_pos = make_pos(2, 2);
    const int initial_bottom_mask = kFullMask;
    const int initial_top_mask = 0;
    return data.expected_nc[static_cast<std::size_t>(initial_bottom_mask)]
                           [static_cast<std::size_t>(initial_top_mask)]
                           [static_cast<std::size_t>(initial_pos)];
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    SolutionData data;
    if (!build_solution_data(data)) {
        std::cerr << "Failed to build linear systems." << '\n';
        return 2;
    }
    if (options.run_checkpoints && !run_checkpoints(data)) {
        return 3;
    }

    const int initial_pos = make_pos(2, 2);
    const ld answer = data.expected_nc[static_cast<std::size_t>(kFullMask)][0]
                                      [static_cast<std::size_t>(initial_pos)];

    std::cout << std::fixed << std::setprecision(6) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

class Grid:
    def __init__(self):
        self.neighbors = [[-1]*4 for _ in range(25)]
        self.degree = [0]*25
        
        for y in range(5):
            for x in range(5):
                p = y * 5 + x
                d = 0
                if x > 0:
                    self.neighbors[p][d] = p - 1
                    d += 1
                if x + 1 < 5:
                    self.neighbors[p][d] = p + 1
                    d += 1
                if y > 0:
                    self.neighbors[p][d] = p - 5
                    d += 1
                if y + 1 < 5:
                    self.neighbors[p][d] = p + 5
                    d += 1
                self.degree[p] = d

def solve_linear_system(a, rhs):
    eps = 1e-15
    for col in range(25):
        pivot = col
        best = abs(a[col][col])
        for row in range(col + 1, 25):
            cand = abs(a[row][col])
            if cand > best:
                best = cand
                pivot = row
                
        if best < eps:
            return False
            
        if pivot != col:
            a[pivot], a[col] = a[col], a[pivot]
            rhs[pivot], rhs[col] = rhs[col], rhs[pivot]
            
        inv_pivot = 1.0 / a[col][col]
        for j in range(col, 25):
            a[col][j] *= inv_pivot
        for r in range(6):
            rhs[col][r] *= inv_pivot
            
        for row in range(25):
            if row == col:
                continue
            factor = a[row][col]
            if abs(factor) < eps:
                continue
            for j in range(col, 25):
                a[row][j] -= factor * a[col][j]
            for r in range(6):
                rhs[row][r] -= factor * rhs[col][r]
                
    return True

class HittingData:
    def __init__(self):
        self.expect = [0.0] * 25
        self.prob = [[0.0]*5 for _ in range(25)]

def solve_hitting_data(grid, row, mask):
    data = HittingData()
    is_target = [False] * 25
    
    for col in range(5):
        if (mask & (1 << col)) != 0:
            is_target[row * 5 + col] = True
            
    a = [[0.0]*25 for _ in range(25)]
    rhs = [[0.0]*6 for _ in range(25)]
    
    for state in range(25):
        if is_target[state]:
            a[state][state] = 1.0
            rhs[state][0] = 0.0
            col = state % 5
            rhs[state][1 + col] = 1.0
            continue
            
        a[state][state] = 1.0
        deg = grid.degree[state]
        step = 1.0 / deg
        for i in range(deg):
            nb = grid.neighbors[state][i]
            a[state][nb] -= step
        rhs[state][0] = 1.0
        
    if not solve_linear_system(a, rhs):
        return None
        
    for state in range(25):
        data.expect[state] = rhs[state][0]
        for col in range(5):
            data.prob[state][col] = rhs[state][1 + col]
            
    return data

def build_solution_data():
    grid = Grid()
    
    hit_top = [None] * 32
    hit_bottom = [None] * 32
    
    for mask in range(1, 32):
        res = solve_hitting_data(grid, 0, mask)
        if res is None:
            return None
        hit_top[mask] = res
        
        res = solve_hitting_data(grid, 4, mask)
        if res is None:
            return None
        hit_bottom[mask] = res
        
    popcount = [bin(m).count('1') for m in range(32)]
    
    expected_nc = [[[0.0]*25 for _ in range(32)] for _ in range(32)]
    expected_c = [[[0.0]*25 for _ in range(32)] for _ in range(32)]
    
    expected_nc[0][31] = [0.0] * 25
    
    for top_count in range(4, -1, -1):
        for top_mask in range(32):
            if popcount[top_mask] != top_count:
                continue
            empty_top_mask = 31 ^ top_mask
            drop = hit_top[empty_top_mask]
            
            for bottom_mask in range(32):
                if popcount[bottom_mask] != 4 - top_count:
                    continue
                    
                for pos in range(25):
                    value = drop.expect[pos]
                    for col in range(5):
                        bit = 1 << col
                        if (empty_top_mask & bit) == 0:
                            continue
                        value += drop.prob[pos][col] * expected_nc[bottom_mask][top_mask | bit][col]
                        
                    expected_c[bottom_mask][top_mask][pos] = value
                    
        for top_mask in range(32):
            if popcount[top_mask] != top_count:
                continue
            for bottom_mask in range(32):
                if popcount[bottom_mask] != 5 - top_count:
                    continue
                    
                pick = hit_bottom[bottom_mask]
                for pos in range(25):
                    value = pick.expect[pos]
                    for col in range(5):
                        bit = 1 << col
                        if (bottom_mask & bit) == 0:
                            continue
                        value += pick.prob[pos][col] * expected_c[bottom_mask ^ bit][top_mask][20 + col]
                        
                    expected_nc[bottom_mask][top_mask][pos] = value
                    
    return expected_nc

def solve():
    expected_nc = build_solution_data()
    if expected_nc is None:
        return "Error"
        
    initial_pos = 2 * 5 + 2
    initial_bottom_mask = 31
    initial_top_mask = 0
    ans = expected_nc[initial_bottom_mask][initial_top_mask][initial_pos]
    return f"{ans:.6f}"

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler280 {
    static final int SIZE = 5;
    static final int CELLS = SIZE * SIZE;
    static final int COLS = 5;
    static final int MASK_COUNT = 1 << COLS;
    static final int FULL_MASK = MASK_COUNT - 1;
    static final int RHS_COUNT = 1 + COLS;

    static class Grid {
        int[][] neighbors = new int[CELLS][4];
        int[] degree = new int[CELLS];

        Grid() {
            for (int y = 0; y < SIZE; ++y) {
                for (int x = 0; x < SIZE; ++x) {
                    int p = y * SIZE + x;
                    int d = 0;
                    if (x > 0)
                        neighbors[p][d++] = p - 1;
                    if (x + 1 < SIZE)
                        neighbors[p][d++] = p + 1;
                    if (y > 0)
                        neighbors[p][d++] = p - SIZE;
                    if (y + 1 < SIZE)
                        neighbors[p][d++] = p + SIZE;
                    degree[p] = d;
                }
            }
        }
    }

    static boolean solveLinearSystem(double[][] a, double[][] rhs) {
        double eps = 1e-15;

        for (int col = 0; col < CELLS; ++col) {
            int pivot = col;
            double best = Math.abs(a[col][col]);
            for (int row = col + 1; row < CELLS; ++row) {
                double cand = Math.abs(a[row][col]);
                if (cand > best) {
                    best = cand;
                    pivot = row;
                }
            }

            if (best < eps)
                return false;

            if (pivot != col) {
                double[] tmpA = a[pivot];
                a[pivot] = a[col];
                a[col] = tmpA;
                double[] tmpR = rhs[pivot];
                rhs[pivot] = rhs[col];
                rhs[col] = tmpR;
            }

            double invPivot = 1.0 / a[col][col];
            for (int j = col; j < CELLS; ++j)
                a[col][j] *= invPivot;
            for (int r = 0; r < RHS_COUNT; ++r)
                rhs[col][r] *= invPivot;

            for (int row = 0; row < CELLS; ++row) {
                if (row == col)
                    continue;
                double factor = a[row][col];
                if (Math.abs(factor) < eps)
                    continue;
                for (int j = col; j < CELLS; ++j)
                    a[row][j] -= factor * a[col][j];
                for (int r = 0; r < RHS_COUNT; ++r)
                    rhs[row][r] -= factor * rhs[col][r];
            }
        }
        return true;
    }

    static class HittingData {
        double[] expect = new double[CELLS];
        double[][] prob = new double[CELLS][COLS];
    }

    static HittingData solveHittingData(Grid grid, int row, int mask) {
        HittingData data = new HittingData();
        boolean[] isTarget = new boolean[CELLS];

        for (int col = 0; col < COLS; ++col) {
            if ((mask & (1 << col)) != 0) {
                isTarget[row * SIZE + col] = true;
            }
        }

        double[][] a = new double[CELLS][CELLS];
        double[][] rhs = new double[CELLS][RHS_COUNT];

        for (int state = 0; state < CELLS; ++state) {
            if (isTarget[state]) {
                a[state][state] = 1.0;
                rhs[state][0] = 0.0;
                int col = state % SIZE;
                rhs[state][1 + col] = 1.0;
                continue;
            }

            a[state][state] = 1.0;
            int deg = grid.degree[state];
            double step = 1.0 / deg;
            for (int i = 0; i < deg; ++i) {
                int nb = grid.neighbors[state][i];
                a[state][nb] -= step;
            }
            rhs[state][0] = 1.0;
        }

        if (!solveLinearSystem(a, rhs))
            return null;

        for (int state = 0; state < CELLS; ++state) {
            data.expect[state] = rhs[state][0];
            for (int col = 0; col < COLS; ++col) {
                data.prob[state][col] = rhs[state][1 + col];
            }
        }

        return data;
    }

    static class SolutionData {
        HittingData[] hitTop = new HittingData[MASK_COUNT];
        HittingData[] hitBottom = new HittingData[MASK_COUNT];
        double[][][] expectedNc = new double[MASK_COUNT][MASK_COUNT][CELLS];
        double[][][] expectedC = new double[MASK_COUNT][MASK_COUNT][CELLS];
    }

    static boolean buildSolutionData(SolutionData out) {
        Grid grid = new Grid();

        for (int mask = 1; mask < MASK_COUNT; ++mask) {
            out.hitTop[mask] = solveHittingData(grid, 0, mask);
            if (out.hitTop[mask] == null)
                return false;
            out.hitBottom[mask] = solveHittingData(grid, SIZE - 1, mask);
            if (out.hitBottom[mask] == null)
                return false;
        }

        int[] popcount = new int[MASK_COUNT];
        for (int mask = 0; mask < MASK_COUNT; ++mask) {
            popcount[mask] = Integer.bitCount(mask);
        }

        for (int pos = 0; pos < CELLS; ++pos) {
            out.expectedNc[0][FULL_MASK][pos] = 0.0;
        }

        for (int topCount = 4; topCount >= 0; --topCount) {
            for (int topMask = 0; topMask < MASK_COUNT; ++topMask) {
                if (popcount[topMask] != topCount)
                    continue;
                int emptyTopMask = FULL_MASK ^ topMask;
                HittingData drop = out.hitTop[emptyTopMask];

                for (int bottomMask = 0; bottomMask < MASK_COUNT; ++bottomMask) {
                    if (popcount[bottomMask] != 4 - topCount)
                        continue;

                    for (int pos = 0; pos < CELLS; ++pos) {
                        double value = drop.expect[pos];
                        for (int col = 0; col < COLS; ++col) {
                            int bit = 1 << col;
                            if ((emptyTopMask & bit) == 0)
                                continue;
                            value += drop.prob[pos][col] * out.expectedNc[bottomMask][topMask | bit][col];
                        }
                        out.expectedC[bottomMask][topMask][pos] = value;
                    }
                }
            }

            for (int topMask = 0; topMask < MASK_COUNT; ++topMask) {
                if (popcount[topMask] != topCount)
                    continue;
                for (int bottomMask = 0; bottomMask < MASK_COUNT; ++bottomMask) {
                    if (popcount[bottomMask] != 5 - topCount)
                        continue;

                    HittingData pick = out.hitBottom[bottomMask];
                    for (int pos = 0; pos < CELLS; ++pos) {
                        double value = pick.expect[pos];
                        for (int col = 0; col < COLS; ++col) {
                            int bit = 1 << col;
                            if ((bottomMask & bit) == 0)
                                continue;
                            value += pick.prob[pos][col]
                                    * out.expectedC[bottomMask ^ bit][topMask][(SIZE - 1) * SIZE + col];
                        }
                        out.expectedNc[bottomMask][topMask][pos] = value;
                    }
                }
            }
        }
        return true;
    }

    public static String solve() {
        SolutionData data = new SolutionData();
        if (!buildSolutionData(data))
            return "Error";

        int initialPos = 2 * SIZE + 2;
        int initialBottomMask = FULL_MASK;
        int initialTopMask = 0;
        double ans = data.expectedNc[initialBottomMask][initialTopMask][initialPos];

        return String.format(Locale.US, "%.6f", ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}