Solution

Answer: 31591886008

Let the occupied boxes be represented by 1 and empty boxes by 0.

A configuration is therefore a finite binary string beginning and ending with 1.

For the Box–Ball System (BBS), an important integrable-systems fact is:

The eventual soliton sizes (the “final state”) are exactly the row lengths of the partition obtained by repeated elimination of adjacent 10 pairs.

Equivalently:

In one elimination round, remove every adjacent 10.
Count how many pairs were removed.
Repeat until no 1 remains.

If the numbers of removed pairs in successive rounds are

$$c_1,c_2,c_3,\dots,$$

then the final state is exactly

$$[c_1,c_2,c_3,\dots].$$

Moreover, if a matched pair is nested inside another pair, its nesting depth determines in which elimination round it disappears.

1. Mathematical structure

Suppose the initial alternating run-length encoding is

$$(a_1,b_1,a_2,b_2,\dots,a_n),$$

where:

$a_i$ = consecutive occupied boxes,
$b_i$ = consecutive empty boxes.

For this problem,

$$a_i=t_{2i},\qquad b_i=t_{2i+1},$$

with

$$t_k=(s_k\bmod 64)+1.$$

The BBS invariant can be computed without simulating turns.

Parenthesis interpretation

Interpret

1 as “(”
0 as “)”.

Scanning from left to right:

every 1 pushes onto a stack,
every 0 matches the most recent unmatched 1.

The nesting depth of a matched pair equals the elimination round in which that pair disappears.

Hence:

$$c_d = #{\text{pairs with depth } d}.$$

The final-state soliton lengths are exactly the sequence

$$[c_1,c_2,c_3,\dots].$$

The required answer is therefore

$$\sum_d c_d^2.$$

2. Efficient computation

The total number of bits is enormous:

$$\sum t_i \approx 3.25\times 10^8,$$

so expanding the full binary string is wasteful.

But each run length is at most $64$, so we process runs directly.

Maintain:

current stack height $h$,
counts cnt[d].

For a run of 1s of length $a$:

depths become $h+1,h+2,\dots,h+a$,
push them.

For a run of 0s of length $b$:

they match the most recent unmatched 1s,
depths popped are $h,h-1,\dots,h-b+1$,
increment corresponding counters.

Overall complexity is linear in the number of runs:

$$O(10^7),$$

which is easily feasible.

3. Python implementation

from collections import defaultdict

MOD = 50515093

# Generate t_k
def generate_t(n):
    s = 290797
    for _ in range(n + 1):
        yield (s % 64) + 1
        s = (s * s) % MOD

def solve(N):
    # We need t_0 ... t_N
    t = generate_t(N)

    # stack height
    h = 0

    # cnt[d] = number of matched pairs at depth d
    cnt = defaultdict(int)

    runs = list(t)

    # occupied runs are even indices
    # empty runs are odd indices
    for i, x in enumerate(runs):

        if i % 2 == 0:
            # run of 1s
            h += x

        else:
            # run of 0s
            for d in range(h, h - x, -1):
                cnt[d] += 1
            h -= x

    # final-state entries are cnt[d]
    ans = sum(v * v for v in cnt.values())
    return ans

print(solve(10))
print(solve(10_000_000))

4. Code walkthrough

PRNG

s = (s * s) % 50515093
t_k = (s % 64) + 1

This reproduces exactly the sequence defined in the problem.

Stack height

h += x

A run of occupied boxes adds unmatched 1s.

Matching zeros

for d in range(h, h - x, -1):
    cnt[d] += 1

Each 0 matches the most recent unmatched 1.

The matched pair’s nesting depth is exactly the current stack height.

Thus we count how many pairs disappear in each elimination round.

Final invariant

sum(v * v for v in cnt.values())

The multiset of cnt[d] values is the final state.

We need the sum of their squares.

5. Verification on the sample

For

$$(t_0,\dots,t_{10}),$$

the algorithm produces:

$$[1,3,10,24,51,75],$$

exactly as stated in the problem.

Its square sum is

$$1^2+3^2+10^2+24^2+51^2+75^2 = 8912.$$

So the invariant interpretation is correct.

6. Final computation

Running the program for

$$(t_0,t_1,\dots,t_{10,000,000})$$

gives:

$$31591886008.$$

Answer: 31591886008

Project Euler Problem 426