Problem 426: Box-Ball System
View on Project EulerProject Euler Problem 426 Solution
EulerSolve provides an optimized solution for Project Euler Problem 426, Box-Ball System, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Set \(s_0=290797\), define $$s_{k+1}=s_k^2 \bmod 50515093,\qquad t_k=(s_k \bmod 64)+1,$$ and build the initial box-ball configuration as the infinite binary word $$W=1^{t_0}0^{t_1}1^{t_2}0^{t_3}\cdots 1^{t_{10^7}}0^\infty,$$ where \(1\) denotes an occupied box and \(0\) an empty box. One BBS turn moves each ball once, from left to right, to the nearest empty box on its right. After sufficiently many turns the configuration separates into stable occupied blocks with lengths \(L_1,\dots,L_m\), and the target quantity is $$\sum_{j=1}^{m} L_j^2.$$ Mathematical Approach Step 1: Replace the Dynamics by Repeated \(10\)-Elimination A standard invariant of the box-ball system is obtained by repeatedly deleting every adjacent pattern \(10\), with all deletions in a given round performed simultaneously and with the infinite tail of zeros already present on the right. Let \(N_r\) be the number of pairs deleted in elimination round \(r\). The stable soliton lengths form the conjugate partition of these row counts. Equivalently, $$N_r=\#\{j:L_j\ge r\}.$$ This identity is the key reduction: once the numbers \(N_r\) are known, the final occupied block lengths are known as well....
Detailed mathematical approach
Problem Summary
Set \(s_0=290797\), define
$$s_{k+1}=s_k^2 \bmod 50515093,\qquad t_k=(s_k \bmod 64)+1,$$
and build the initial box-ball configuration as the infinite binary word
$$W=1^{t_0}0^{t_1}1^{t_2}0^{t_3}\cdots 1^{t_{10^7}}0^\infty,$$
where \(1\) denotes an occupied box and \(0\) an empty box. One BBS turn moves each ball once, from left to right, to the nearest empty box on its right. After sufficiently many turns the configuration separates into stable occupied blocks with lengths \(L_1,\dots,L_m\), and the target quantity is
$$\sum_{j=1}^{m} L_j^2.$$
Mathematical Approach
Step 1: Replace the Dynamics by Repeated \(10\)-Elimination
A standard invariant of the box-ball system is obtained by repeatedly deleting every adjacent pattern \(10\), with all deletions in a given round performed simultaneously and with the infinite tail of zeros already present on the right. Let \(N_r\) be the number of pairs deleted in elimination round \(r\).
The stable soliton lengths form the conjugate partition of these row counts. Equivalently,
$$N_r=\#\{j:L_j\ge r\}.$$
This identity is the key reduction: once the numbers \(N_r\) are known, the final occupied block lengths are known as well. The number of blocks of exact length \(r\) is therefore
$$N_r-N_{r+1}.$$
Geometrically, each final block of length \(\ell\) contributes one cell to each of the first \(\ell\) elimination rows, so the Ferrers diagram of the stable blocks has row lengths \(N_1,N_2,\dots\).
Step 2: A One-Pass Stack Interpretation
The implementations compute the elimination rounds without physically deleting symbols. Scan the binary word from left to right and keep a stack of unmatched occupied boxes. For each unmatched \(1\), store the largest elimination round of any matched pair nested inside it.
When an occupied box is read, push \(0\): at that moment no inner pair has been formed yet. When an empty box is read, it matches the most recent unmatched occupied box. If the stored maximum for that occupied box is \(d\), then all nested pairs vanish by round \(d\), so the outer pair becomes adjacent exactly one round later. Its elimination round is therefore
$$d+1.$$
After recording one more pair in round \(d+1\), this value must influence the nearest unmatched occupied box on the left, because that box now contains a nested structure whose maximum round may have increased. Thus the new value is propagated leftward by a max-update on the stack top.
This is the same recurrence as parenthesis matching with nesting depth: a pair disappears one round after everything strictly inside it has disappeared. At the end of the explicit word, the infinite zero tail is handled by draining the remaining stack in the same way.
Step 3: Work Directly with the Run-Length Encoding
The initial state is given as alternating run lengths, so the algorithm never materializes a full turn-by-turn BBS trajectory. For an occupied run of length \(\ell\), it pushes \(\ell\) zeros onto the stack. For an empty run of length \(\ell\), it performs up to \(\ell\) closures, stopping early if no unmatched occupied box remains.
Because every \(t_k\) lies in \(\{1,\dots,64\}\), each generated run costs only \(O(1)\) work, and the complete computation is just a streaming pass over the pseudorandom sequence.
Step 4: Convert Row Counts into the Required Sum
Once the counts \(N_r\) are available, the sum of squares follows from the elementary identity
$$\ell^2=\sum_{r=1}^{\ell}(2r-1).$$
Summing over all stable blocks and exchanging the order of summation gives
$$\sum_{j=1}^{m}L_j^2=\sum_{j=1}^{m}\sum_{r=1}^{L_j}(2r-1)=\sum_{r\ge 1}(2r-1)\#\{j:L_j\ge r\}.$$
Now substitute \(N_r=\#\{j:L_j\ge r\}\) to obtain the exact formula used by the implementations:
$$\boxed{\sum_{j=1}^{m}L_j^2=\sum_{r\ge 1}(2r-1)N_r.}$$
Worked Examples
Consider the small run encoding \([2,2,2,1,2]\). It represents the infinite word
$$110011011\,0^\infty.$$
Repeated \(10\)-elimination removes \(3\) pairs in the first round, \(2\) in the second, and \(1\) in the third. Hence
$$N_1=3,\qquad N_2=2,\qquad N_3=1.$$
The exact multiplicities are
$$N_1-N_2=1,\qquad N_2-N_3=1,\qquad N_3-N_4=1,$$
so the stable occupied blocks have lengths \([1,2,3]\).
For the \(11\)-run sample from the problem statement, the stable lengths are \([1,3,10,24,51,75]\). Therefore \(N_r\) is piecewise constant:
$$N_r=\begin{cases} 6,& r=1,\\ 5,& 2\le r\le 3,\\ 4,& 4\le r\le 10,\\ 3,& 11\le r\le 24,\\ 2,& 25\le r\le 51,\\ 1,& 52\le r\le 75,\\ 0,& r\ge 76. \end{cases}$$
Substituting these row counts into the boxed identity yields
$$\sum_{r\ge 1}(2r-1)N_r=8912,$$
which matches the sample value exactly.
How the Code Works
The C++, Python, and Java implementations all stream the pseudorandom generator, alternate between occupied and empty runs, and maintain the same stack invariant. Every closure contributes one unit to the appropriate elimination round, and a final drain accounts for the infinite empty tail. After that single pass, the implementations evaluate \(\sum_{r\ge 1}(2r-1)N_r\) directly from the stored round counts. No explicit simulation of many BBS turns is performed.
Complexity Analysis
Let \(B\) be the total number of occupied boxes in the encoded initial state. Each occupied box is pushed onto the stack once and popped once, so the total stack work is \(O(B)\). The round-count array has length equal to the largest elimination round encountered. Thus the memory usage is \(O(H+D)\), where \(H\) is the maximum active stack height and \(D\) the maximum round index; both are bounded by \(B\).
Because every run length satisfies \(1\le t_k\le 64\), we also have \(B=O(R)\) for \(R=10^7+1\) generated runs. Therefore the full algorithm is linear in the size of the encoded input and dramatically faster than any direct simulation over many BBS turns.
Footnotes and References
- Problem page: https://projecteuler.net/problem=426
- D. Takahashi and J. Satsuma, A Soliton Cellular Automaton, Journal of the Physical Society of Japan, 59(10), 3514-3519, 1990.
- Box-ball system overview: Wikipedia — Box-ball system
- Partition conjugation and Ferrers diagrams: Wikipedia — Partition (number theory)
- Run-length encoding: Wikipedia — Run-length encoding
Problem 426 source code
C++
#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <mutex>
#include <string>
#include <thread>
#include <unordered_set>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr int kTargetIndex = 10'000'000;
constexpr u64 kSeed = 290'797ULL;
constexpr u64 kMod = 50'515'093ULL;
std::string to_string_u128(u128 value) {
if (value == 0) {
return "0";
}
std::string digits;
while (value > 0) {
const int digit = static_cast<int>(value % 10);
digits.push_back(static_cast<char>('0' + digit));
value /= 10;
}
std::reverse(digits.begin(), digits.end());
return digits;
}
std::vector<int> generate_t_runs(int max_index) {
std::vector<int> runs;
runs.reserve(static_cast<std::size_t>(max_index) + 1);
u64 s = kSeed;
for (int i = 0; i <= max_index; ++i) {
runs.push_back(static_cast<int>(s % 64ULL) + 1);
s = static_cast<u64>((static_cast<u128>(s) * s) % kMod);
}
return runs;
}
class InvariantProfileSolver {
public:
void process_run(int run_length, bool occupied) {
if (run_length <= 0) {
return;
}
if (occupied) {
stack_.resize(stack_.size() + static_cast<std::size_t>(run_length), 0U);
return;
}
while (run_length > 0 && !stack_.empty()) {
close_top_pair();
--run_length;
}
}
void process_runs(const std::vector<int>& runs) {
reset();
bool occupied = true;
for (const int run_length : runs) {
process_run(run_length, occupied);
occupied = !occupied;
}
finalize();
}
void finalize() {
while (!stack_.empty()) {
close_top_pair();
}
}
u128 sum_of_square_final_state() const {
u128 answer = 0;
for (std::size_t round = 1; round < round_counts_.size(); ++round) {
answer += static_cast<u128>(2 * round - 1) * static_cast<u128>(round_counts_[round]);
}
return answer;
}
std::vector<u64> final_state_lengths() const {
// This is used only in small validation checkpoints.
std::vector<u64> lengths;
for (std::size_t len = 1; len < round_counts_.size(); ++len) {
const u64 current = round_counts_[len];
const u64 next = (len + 1 < round_counts_.size()) ? round_counts_[len + 1] : 0;
if (current < next) {
return {};
}
const u64 multiplicity = current - next;
lengths.insert(lengths.end(), static_cast<std::size_t>(multiplicity),
static_cast<u64>(len));
}
return lengths;
}
private:
std::vector<u32> stack_;
std::vector<u64> round_counts_{0}; // 1-indexed by elimination round.
void reset() {
stack_.clear();
round_counts_.assign(1, 0);
}
void close_top_pair() {
const u32 child_max_round = stack_.back();
stack_.pop_back();
const std::size_t round = static_cast<std::size_t>(child_max_round) + 1;
if (round >= round_counts_.size()) {
round_counts_.resize(round + 1, 0);
}
++round_counts_[round];
if (!stack_.empty() && stack_.back() < round) {
stack_.back() = static_cast<u32>(round);
}
}
};
std::vector<int> runs_to_positions(const std::vector<int>& runs) {
std::vector<int> positions;
int cursor = 0;
bool occupied = true;
for (const int run_length : runs) {
if (occupied) {
positions.reserve(positions.size() + static_cast<std::size_t>(run_length));
for (int offset = 0; offset < run_length; ++offset) {
positions.push_back(cursor + offset);
}
}
cursor += run_length;
occupied = !occupied;
}
return positions;
}
std::vector<int> positions_to_runs(const std::vector<int>& positions) {
if (positions.empty()) {
return {};
}
std::vector<int> runs;
int i = 0;
while (i < static_cast<int>(positions.size())) {
int j = i + 1;
while (j < static_cast<int>(positions.size()) &&
positions[static_cast<std::size_t>(j)] ==
positions[static_cast<std::size_t>(j - 1)] + 1) {
++j;
}
runs.push_back(j - i);
if (j < static_cast<int>(positions.size())) {
const int gap = positions[static_cast<std::size_t>(j)] -
positions[static_cast<std::size_t>(j - 1)] - 1;
runs.push_back(gap);
}
i = j;
}
return runs;
}
void apply_one_turn(std::vector<int>& positions) {
if (positions.empty()) {
return;
}
std::unordered_set<int> occupied;
occupied.reserve(positions.size() * 4 + 8);
for (const int position : positions) {
occupied.insert(position);
}
for (std::size_t i = 0; i < positions.size(); ++i) {
const int current = positions[i];
occupied.erase(current);
int target = current + 1;
while (occupied.find(target) != occupied.end()) {
++target;
}
occupied.insert(target);
positions[i] = target;
}
}
std::vector<int> occupied_block_lengths(const std::vector<int>& positions) {
if (positions.empty()) {
return {};
}
std::vector<int> blocks;
int i = 0;
while (i < static_cast<int>(positions.size())) {
int j = i + 1;
while (j < static_cast<int>(positions.size()) &&
positions[static_cast<std::size_t>(j)] ==
positions[static_cast<std::size_t>(j - 1)] + 1) {
++j;
}
blocks.push_back(j - i);
i = j;
}
return blocks;
}
std::vector<int> brute_force_final_blocks_after_turns(const std::vector<int>& runs,
int turns) {
std::vector<int> positions = runs_to_positions(runs);
for (int t = 0; t < turns; ++t) {
apply_one_turn(positions);
}
return occupied_block_lengths(positions);
}
std::vector<u64> fast_final_blocks(const std::vector<int>& runs) {
InvariantProfileSolver solver;
solver.process_runs(runs);
return solver.final_state_lengths();
}
template <typename T>
std::string vector_to_string(const std::vector<T>& values) {
std::string out = "[";
for (std::size_t i = 0; i < values.size(); ++i) {
if (i > 0) {
out += ", ";
}
out += std::to_string(values[i]);
}
out += "]";
return out;
}
bool run_turn_rule_checkpoint() {
const std::vector<int> initial = {2, 2, 2, 1, 2};
const std::vector<int> expected = {2, 2, 1, 2, 3};
std::vector<int> positions = runs_to_positions(initial);
apply_one_turn(positions);
const std::vector<int> got = positions_to_runs(positions);
if (got != expected) {
std::cerr << "Turn checkpoint failed: expected " << vector_to_string(expected)
<< ", got " << vector_to_string(got) << '\n';
return false;
}
return true;
}
bool run_statement_sample_checkpoint() {
const std::vector<int> sample_runs = generate_t_runs(10);
InvariantProfileSolver solver;
solver.process_runs(sample_runs);
const std::vector<u64> expected_lengths = {1, 3, 10, 24, 51, 75};
const std::vector<u64> got_lengths = solver.final_state_lengths();
if (got_lengths != expected_lengths) {
std::cerr << "Sample checkpoint failed: expected final state "
<< vector_to_string(expected_lengths) << ", got "
<< vector_to_string(got_lengths) << '\n';
return false;
}
const u64 expected_sum = 8'912;
const u64 got_sum = static_cast<u64>(solver.sum_of_square_final_state());
if (got_sum != expected_sum) {
std::cerr << "Sample sum-of-squares checkpoint failed: expected " << expected_sum
<< ", got " << got_sum << '\n';
return false;
}
return true;
}
struct SmallCase {
std::vector<int> runs;
};
std::vector<SmallCase> build_small_cross_check_cases() {
constexpr int kCaseCount = 96;
std::vector<SmallCase> cases;
cases.reserve(kCaseCount);
u64 state = 0x9e3779b97f4a7c15ULL;
auto next_u32 = [&state]() {
state = state * 6364136223846793005ULL + 1442695040888963407ULL;
return static_cast<u32>(state >> 32);
};
for (int cid = 0; cid < kCaseCount; ++cid) {
const int occupied_runs = 1 + static_cast<int>(next_u32() % 6U);
SmallCase test_case;
test_case.runs.reserve(static_cast<std::size_t>(2 * occupied_runs - 1));
for (int i = 0; i < 2 * occupied_runs - 1; ++i) {
test_case.runs.push_back(1 + static_cast<int>(next_u32() % 6U));
}
cases.push_back(std::move(test_case));
}
return cases;
}
bool run_small_cross_checks(unsigned requested_threads) {
constexpr int kBruteTurns = 2'000;
const std::vector<SmallCase> cases = build_small_cross_check_cases();
unsigned threads = requested_threads;
if (threads == 0) {
threads = 1;
}
threads = std::max(1u, std::min(threads, static_cast<unsigned>(cases.size())));
std::atomic<std::size_t> next_case{0};
std::atomic<bool> failed{false};
std::mutex failure_mutex;
std::string failure_message;
auto worker = [&]() {
while (!failed.load(std::memory_order_relaxed)) {
const std::size_t index = next_case.fetch_add(1, std::memory_order_relaxed);
if (index >= cases.size()) {
break;
}
const std::vector<u64> fast_u64 = fast_final_blocks(cases[index].runs);
const std::vector<int> slow =
brute_force_final_blocks_after_turns(cases[index].runs, kBruteTurns);
if (fast_u64.size() != slow.size()) {
if (!failed.exchange(true, std::memory_order_relaxed)) {
std::lock_guard<std::mutex> lock(failure_mutex);
failure_message =
"Random cross-check failed on case " + std::to_string(index) +
": expected " + vector_to_string(slow) + ", got " +
vector_to_string(fast_u64) + ", runs=" +
vector_to_string(cases[index].runs);
}
break;
}
bool mismatch = false;
for (std::size_t i = 0; i < slow.size(); ++i) {
if (static_cast<u64>(slow[i]) != fast_u64[i]) {
mismatch = true;
break;
}
}
if (mismatch) {
if (!failed.exchange(true, std::memory_order_relaxed)) {
std::lock_guard<std::mutex> lock(failure_mutex);
failure_message =
"Random cross-check failed on case " + std::to_string(index) +
": expected " + vector_to_string(slow) + ", got " +
vector_to_string(fast_u64) + ", runs=" +
vector_to_string(cases[index].runs);
}
break;
}
}
};
if (threads == 1) {
worker();
} else {
std::vector<std::thread> pool;
pool.reserve(threads);
for (unsigned t = 0; t < threads; ++t) {
pool.emplace_back(worker);
}
for (std::thread& thread : pool) {
thread.join();
}
}
if (failed.load(std::memory_order_relaxed)) {
std::cerr << failure_message << '\n';
return false;
}
return true;
}
bool run_validation_checkpoints(unsigned threads) {
if (!run_turn_rule_checkpoint()) {
return false;
}
if (!run_statement_sample_checkpoint()) {
return false;
}
if (!run_small_cross_checks(threads)) {
return false;
}
return true;
}
u128 solve_target() {
InvariantProfileSolver solver;
u64 s = kSeed;
bool occupied = true;
for (int index = 0; index <= kTargetIndex; ++index) {
const int run_length = static_cast<int>(s % 64ULL) + 1;
solver.process_run(run_length, occupied);
occupied = !occupied;
s = static_cast<u64>((static_cast<u128>(s) * s) % kMod);
}
solver.finalize();
return solver.sum_of_square_final_state();
}
} // namespace
int main() {
unsigned threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 4;
}
if (!run_validation_checkpoints(threads)) {
return 1;
}
const u128 answer = solve_target();
std::cout << to_string_u128(answer) << '\n';
return 0;
}
Python
def solve():
SEED = 290797; MOD = 50515093; TARGET = 10000000
class Solver:
def __init__(self):
self.stack = []; self.rc = [0]
def process_run(self, rl, occ):
if rl <= 0: return
if occ: self.stack.extend([0]*rl); return
for _ in range(rl):
if not self.stack: break
cmr = self.stack.pop()
rd = cmr + 1
while len(self.rc) <= rd: self.rc.append(0)
self.rc[rd] += 1
if self.stack and self.stack[-1] < rd: self.stack[-1] = rd
def finalize(self):
while self.stack:
cmr = self.stack.pop()
rd = cmr + 1
while len(self.rc) <= rd: self.rc.append(0)
self.rc[rd] += 1
if self.stack and self.stack[-1] < rd: self.stack[-1] = rd
def sos(self):
ans = 0
for r in range(1, len(self.rc)):
ans += (2*r-1)*self.rc[r]
return ans
sol = Solver()
s = SEED; occ = True
for _ in range(TARGET+1):
rl = int(s % 64) + 1
sol.process_run(rl, occ)
occ = not occ
s = (s*s) % MOD
sol.finalize()
return str(sol.sos())
if __name__ == '__main__':
print(solve())
Java
public class Euler426 {
private static final int TARGET_INDEX = 10000000;
private static final long SEED = 290797L;
private static final long MOD = 50515093L;
public static String solve() {
int[] stack = new int[TARGET_INDEX + 100];
int stackSize = 0;
long[] roundCounts = new long[TARGET_INDEX / 2 + 100];
int maxRound = 0;
long s = SEED;
boolean occupied = true;
for (int i = 0; i <= TARGET_INDEX; i++) {
int runLength = (int) (s % 64) + 1;
if (occupied) {
while (stackSize + runLength > stack.length) {
int[] newStack = new int[stack.length * 2];
System.arraycopy(stack, 0, newStack, 0, stackSize);
stack = newStack;
}
for (int j = 0; j < runLength; j++) {
stack[stackSize++] = 0;
}
} else {
int rl = runLength;
while (rl > 0 && stackSize > 0) {
int childMaxRound = stack[--stackSize];
int rnd = childMaxRound + 1;
if (rnd > maxRound) {
maxRound = rnd;
if (maxRound >= roundCounts.length) {
long[] newCounts = new long[Math.max(roundCounts.length * 2, maxRound + 100)];
System.arraycopy(roundCounts, 0, newCounts, 0, roundCounts.length);
roundCounts = newCounts;
}
}
roundCounts[rnd]++;
if (stackSize > 0 && stack[stackSize - 1] < rnd) {
stack[stackSize - 1] = rnd;
}
rl--;
}
}
occupied = !occupied;
s = (s * s) % MOD;
}
while (stackSize > 0) {
int childMaxRound = stack[--stackSize];
int rnd = childMaxRound + 1;
if (rnd > maxRound) {
maxRound = rnd;
if (maxRound >= roundCounts.length) {
long[] newCounts = new long[Math.max(roundCounts.length * 2, maxRound + 100)];
System.arraycopy(roundCounts, 0, newCounts, 0, roundCounts.length);
roundCounts = newCounts;
}
}
roundCounts[rnd]++;
if (stackSize > 0 && stack[stackSize - 1] < rnd) {
stack[stackSize - 1] = rnd;
}
}
long ans = 0;
for (int rnd = 1; rnd <= maxRound; rnd++) {
ans += (2L * rnd - 1) * roundCounts[rnd];
}
return String.valueOf(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}