Problem 280: Ant and Seeds
View on Project EulerProject 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 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.
2) The popcount invariant, and how many valid compressed states exist
When the ant is not carrying, every undelivered seed is still sitting at the bottom, and every delivered seed already occupies one top slot. Therefore
$$|B|+|T|=5.$$
When the ant is carrying, one seed has already been removed from the bottom but has not yet been placed on the top, so
$$|B|+|T|=4.$$
This gives the exact set of valid compressed states. Their counts are
$$25\sum_{t=0}^{5}\binom{5}{t}\binom{5}{5-t}=25\binom{10}{5}=6300$$
for non-carrying states, and
$$25\sum_{t=0}^{4}\binom{5}{t}\binom{5}{4-t}=25\binom{10}{4}=5250$$
for carrying states, for a total of
$$6300+5250=11550$$
compressed Markov states. The code does not solve one giant \(11550\times11550\) linear system. It uses a better decomposition.
3) The key idea: jump from one meaningful event to the next
Suppose the current state is \((B,T,p,\text{not carrying})\). Then the next meaningful event is:
the first time the random walk hits one of the active bottom cells
$$bot(B)=\{b_j:\ j\in B\},$$
where \(b_j\) is the bottom cell in column \(j\).
Similarly, in a carrying state \((B,T,p,\text{carrying})\), the next meaningful event is the first time the walk hits one of the still-empty top cells
$$top(\overline T)=\{t_j:\ j\notin T\},$$
where \(t_j\) is the top cell in column \(j\).
This jump is exact, not heuristic. Let \(\tau\) be that first hitting time. By the strong Markov property of the random walk, once the process reaches the event cell, the future depends only on the new compressed state, not on the detailed path used to get there. So we may split the expectation into:
expected travel time to the next event
plus
expected continuation value from the event state.
4) Hitting-time and first-hit-column data
Fix any target set \(S\) contained in the top row or bottom row. For each start cell \(u\), define
$$H_u^{S}=\mathbb E_u[\tau_S],$$
where \(\tau_S\) is the first time the walk hits \(S\). Also define, for each target column \(j\),
$$P_{u\to j}^{S}=\Pr_u(\text{the first hit in }S\text{ occurs at the target cell in column }j).$$
These are determined by first-step decomposition. If \(u\in S\), then the hitting time is already zero. If \(u\notin S\), the first move consumes one step and then averages over neighbors:
$$H_u^{S}=\begin{cases} 0,&u\in S,\\ 1+\dfrac1{\deg(u)}\sum_{v\sim u}H_v^{S},&u\notin S. \end{cases}$$
Likewise, the first-hit-column probabilities satisfy
$$P_{u\to j}^{S}=\begin{cases} 1,&u\text{ is the target cell of column }j,\\ 0,&u\in S\text{ but in a different target column},\\ \dfrac1{\deg(u)}\sum_{v\sim u}P_{v\to j}^{S},&u\notin S. \end{cases}$$
The important structural point is that the coefficient matrix is the same for all right-hand sides. So for each mask the code solves one \(25\times25\) system with 6 right-hand sides: one for the expected time and five for the first-hit probabilities of the five columns.
5) Event-level dynamic programming recurrences
Let \(E_{nc}(B,T,p)\) be the expected remaining number of steps from a non-carrying state, and \(E_c(B,T,p)\) the corresponding quantity from a carrying state.
If the ant is not carrying, the next event is the first hit of \(bot(B)\). Therefore
$$E_{nc}(B,T,p)=H_p^{bot(B)}+\sum_{j\in B}P_{p\to j}^{bot(B)}\,E_c(B\setminus\{j\},T,b_j).$$
The term \(H_p^{bot(B)}\) counts the expected number of ordinary walk steps until pickup. The sum averages over which bottom seed is picked up first. After reaching \(b_j\), the seed is removed from the bottom immediately, so the next state is carrying with bottom mask \(B\setminus\{j\}\).
If the ant is carrying, the next event is the first hit of an empty top slot. Hence
$$E_c(B,T,p)=H_p^{top(\overline T)}+\sum_{j\notin T}P_{p\to j}^{top(\overline T)}\,E_{nc}(B,T\cup\{j\},t_j).$$
Again, the first term is the expected travel time until drop, and the sum averages over which empty top column is reached first. Once \(t_j\) is reached, that top slot becomes occupied immediately.
6) Why this DP has a clean evaluation order
The carrying recurrence increases \(|T|\) by one:
$$E_c(\cdots,T,\cdots)\longrightarrow E_{nc}(\cdots,T\cup\{j\},\cdots).$$
The non-carrying recurrence keeps \(|T|\) unchanged:
$$E_{nc}(\cdots,T,\cdots)\longrightarrow E_c(\cdots,T,\cdots).$$
So if we process states by
$$top\_count=|T|$$
from \(4\) down to \(0\), then:
carrying states at level \(top\_count\) depend only on non-carrying states at level \(top\_count+1\), which are already known;
non-carrying states at level \(top\_count\) depend only on carrying states at the same level, which have just been computed.
That is why the code has a simple layer-by-layer order and never needs iterative relaxation.
7) Base case and initial state
When all five seeds have already been delivered, we must be in a non-carrying state with
$$B=\varnothing,\qquad T=\{0,1,2,3,4\}.$$
No further steps are needed, so
$$E_{nc}(\varnothing,\text{full},p)=0$$
for every position \(p\).
The initial state is
$$B=\text{full},\qquad T=0,\qquad p=\text{center},\qquad \text{not carrying},$$
so the answer is
$$E_{nc}(\text{full},0,\text{center}).$$
8) Worked one-seed sanity check
Suppose only one bottom seed remains, in column \(s\), and only one top slot is still empty, in column \(t\). Then there is no branching left.
If the ant is non-carrying at position \(p\), it must first hit \(b_s\), pick up the seed, and then hit \(t_t\). Therefore
$$E_{nc}(\{s\},\text{full}\setminus\{t\},p)=H_p^{\{b_s\}}+H_{b_s}^{\{t_t\}}.$$
If the ant is already carrying, then only the second part remains:
$$E_c(\varnothing,\text{full}\setminus\{t\},p)=H_p^{\{t_t\}}.$$
The program checks exactly these identities numerically for every choice of \(s\), \(t\), and \(p\). This is a strong checkpoint because it tests both the hitting systems and the event-level DP.
Code Logic
1) Grid construction. build_grid() stores the neighbors and degrees of the 25 cells.
2) Linear algebra core. solve_hitting_data() builds the common coefficient matrix for a target mask on the top or bottom row. solve_linear_system() performs Gaussian elimination with partial pivoting.
3) Precomputation. For every nonempty mask, the code stores hitting-time and first-hit-column data for top targets and bottom targets separately.
4) DP tables. expected_nc[bottom_mask][top_mask][pos] and expected_c[...] implement the two recurrences above.
5) Layer order. The loop over top_count runs from 4 down to 0. For each layer the carrying table is filled first, then the non-carrying table.
6) Checkpoints. The code verifies two invariants: the hit probabilities for an active target mask always sum to \(1\), and the one-seed formulas above agree with the DP tables.
Complexity Analysis
There are only \(31\) nonempty masks per row that need hitting data. Each one solves a \(25\times25\) linear system with 6 right-hand sides, which is tiny. After that, the dynamic program visits only the valid \((B,T,p)\) states described above. So the computation is exact enough for floating-point arithmetic and dramatically smaller than solving a single giant Markov-chain system on the full state space.
Further Reading
- Problem page: https://projecteuler.net/problem=280
- Absorbing Markov chains: https://en.wikipedia.org/wiki/Absorbing_Markov_chain
- Random-walk hitting times: https://en.wikipedia.org/wiki/Random_walk
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());
}
}