Detailed mathematical approach

Problem Summary

For each pair \(1 \le c,d \le 160\), let \(s=c+d\). We place nonnegative integers on a cycle of length \(s\), start from a single pulse \(k\) at one position, and then repeatedly perform a cyclic in-place averaging update with a floor. The key quantity is \(K(c,d)\), the largest starting value that still eventually dies out to the all-zero state.

Once \(K(c,d)\) is known, the problem asks for

\[ G(c,d)=2K(c,d)+1, \qquad \sum_{c=1}^{160}\sum_{d=1}^{160} G(c,d). \]

A brute-force search over all starting values is impossible because the thresholds become large. The implementations succeed by exploiting the monotone structure of the update rule, a reduction by \(\gcd(c,d)\), and a further collapse of the reduced problem to one-dimensional threshold tables.

Mathematical Approach

The whole problem is governed by a discrete dynamical system on a cycle. The useful mathematics is not an abstract template; it is the exact recurrence, the exact invariant regions, and the exact threshold structure used by the search.

The cyclic averaging recurrence

Index the cycle by \(0,1,\dots,s-1\), and let the state after some number of elementary updates be \(x=(x_0,\dots,x_{s-1})\). The initial state is

\[ x_{s-1}=k, \qquad x_i=0 \text{ for } i\ne s-1. \]

At each elementary step, only one coordinate changes. If the current head position is \(h\), then the update is

\[ x_h \leftarrow \left\lfloor \frac{x_h + x_{h-d \bmod s}}{2} \right\rfloor, \]

and then the head advances to \(h+1 \bmod s\). Because \(s=c+d\), the same rule can also be written as

\[ x_h \leftarrow \left\lfloor \frac{x_h + x_{h+c \bmod s}}{2} \right\rfloor. \]

This is the exact recurrence implemented in all three languages. The process is asynchronous: values written earlier in the sweep are immediately visible to later updates in the same sweep.

Order preservation and absorbing regions

Two facts drive the search.

First, the local map \((a,b)\mapsto \lfloor(a+b)/2\rfloor\) is monotone in both inputs. So if two initial pulses satisfy \(k_1 \le k_2\), and we run the same update schedule on both, then the componentwise order is preserved forever. Larger initial data can never produce a smaller trajectory than smaller initial data.

Second, the dynamics has two obvious forward-invariant regions:

\[ x\equiv 0 \implies \text{all future states are } 0, \]

and

\[ x_i \ge 1 \text{ for every } i \implies \left\lfloor \frac{x_i+x_j}{2} \right\rfloor \ge 1, \]

so a state with no zeros can never reach the all-zero state later. This is why the implementations can stop as soon as either every entry is zero or every entry is positive.

There is also a simple boundedness invariant:

\[ 0 \le x_i \le k \qquad \text{for every coordinate and every time.} \]

No update can create a negative value or exceed the initial pulse.

The extinction threshold \(K(c,d)\)

Define the extinction predicate

\[ \mathcal{E}_{c,d}(k)= \begin{cases} 1, & \text{if the trajectory started from } k \text{ reaches } 0,\\ 0, & \text{otherwise.} \end{cases} \]

By monotonicity, \(\mathcal{E}_{c,d}(k)\) is non-increasing in \(k\). Therefore the set of extinguishing initial values is an initial interval

\[ \{k\ge 0 : \mathcal{E}_{c,d}(k)=1\}=\{0,1,\dots,K(c,d)\}. \]

This is the mathematical reason binary search is valid: once some \(k\) fails to die out, every larger \(k\) also fails.

Reduction by the greatest common divisor

Let

\[ g=\gcd(c,d), \qquad r_c=\frac{c}{g}, \qquad r_d=\frac{d}{g}. \]

Since the update for position \(h\) reads from \(h+c \bmod s\), every dependency preserves the residue class modulo \(g\). So the cycle of length \(s=c+d\) splits into \(g\) independent smaller cycles, each of length

\[ \frac{s}{g}=r_c+r_d. \]

Only one of those residue classes contains the initial pulse. The other \(g-1\) classes start at zero and stay zero forever. Therefore the extinction threshold depends only on the reduced pair:

\[ K(c,d)=K(r_c,r_d). \]

This reduction is exact, not heuristic. It explains why the final summation always begins by dividing \(c\) and \(d\) by their gcd.

The one-parameter collapse used by the solver

After the gcd reduction, the implementations do not build a full two-dimensional table of reduced thresholds. Instead they use a stronger structural relation and reconstruct every reduced pair from one-dimensional families.

Write

\[ K_1(n)=K(1,n), \qquad K_{d=1}(m)=K(m,1). \]

The reduced threshold used by the solver is then

\[ K(r_c,r_d)= \begin{cases} K_{d=1}(r_c), & r_d=1,\\ K_1\!\left(r_d+\left\lfloor\frac{r_c-1}{2}\right\rfloor\right), & r_d>1. \end{cases} \]

So the two-parameter family is collapsed to the line \(K(1,n)\), plus a small exceptional strip for \(K(m,1)\). In particular, once \(K_1(1),K_1(2),\dots,K_1(239)\) are known, almost every reduced pair needed by the problem is already determined.

Worked example: the threshold for \((c,d)=(2,1)\)

Here \(s=3\), so the initial state is \((0,0,k)\). The update repeatedly averages a position with the previous one on the 3-cycle.

For \(k=3\), one possible compressed trace is

\[ (0,0,3)\to(1,0,3)\to(1,0,1)\to(1,0,0)\to(0,0,0). \]

So \(k=3\) does become extinct.

For \(k=4\), the process behaves differently:

\[ (0,0,4)\to(2,0,4)\to(2,1,4)\to(2,1,2)\to(2,1,1)\to(1,1,1). \]

Once the state is \((1,1,1)\), every coordinate is positive, so extinction is impossible from then on. Hence \(K(2,1)=3\), and therefore

\[ G(2,1)=2\cdot 3+1=7. \]

This small example captures the general pattern: monotonicity turns the dynamical question into a threshold search.

How the Code Works

Exact extinction tests

The C++, Python, and Java implementations simulate the in-place recurrence exactly for a fixed triple \((c,d,k)\). They maintain the cyclic state, the current head, and a counter of how many entries are zero. That counter gives two exact stopping conditions: all zero means extinction has been reached, while no zero means the state has entered the forward-invariant positive region and cannot die out anymore.

Precomputing reduced thresholds

Each exact threshold is found by doubling an upper bound until extinction fails, then binary-searching between the last extinguishing value and the first non-extinguishing one. For the special family \(K_1(n)=K(1,n)\), the implementations use an eight-way cubic predictor indexed by \(n \bmod 8\) to place the initial search window near the true answer. That predictor changes only the speed of the search, not its correctness.

The one-dimensional table \(K_1(1),\dots,K_1(239)\) is the main precomputation. The C++ and Java implementations parallelize it across threads; the Python implementation uses multiple processes when available and falls back to serial evaluation if needed. For the branch \(K(m,1)\), a few small values are computed directly, and the larger ones are reconstructed from the same collapsed family with direct spot checks for safety.

Final reconstruction of the sum

For every pair \((c,d)\), the implementation first reduces it to \((r_c,r_d)\) by dividing by \(\gcd(c,d)\). It then recovers the needed threshold from the precomputed one-dimensional tables, forms

\[ G(c,d)=2K(r_c,r_d)+1, \]

and accumulates the result over all \(160^2\) pairs. Because the final answer is far larger than a machine word, the total is stored with arbitrary-precision integer arithmetic.

Complexity Analysis

Let \(\tau(c,d,k)\) be the number of elementary updates performed before the extinction test for \((c,d,k)\) reaches one of its two exact stopping conditions. Then one threshold query costs

\[ O\!\big(\tau(c,d,\cdot)\log K(c,d)\big), \]

because doubling plus binary search makes only logarithmically many extinction tests.

The dominant work is the precomputation of the table \(K_1(1),\dots,K_1(239)\). After that, the final double sum over all \(c,d\le 160\) is cheap: each pair needs one gcd, one table lookup, and a few arithmetic operations. Memory usage is small. A single simulation stores at most \(c+d\le 320\) integers, and the precomputed tables themselves have only a few hundred entries.

Footnotes and References

1. Problem page: Project Euler 986
2. Greatest common divisor: Wikipedia - Greatest common divisor
3. Binary search algorithm: Wikipedia - Binary search algorithm
4. Discrete dynamical system: Wikipedia - Discrete dynamical system
5. Asynchronous iteration: Wikipedia - Asynchronous iteration
6. Floor and ceiling functions: Wikipedia - Floor and ceiling functions