Detailed mathematical approach

Problem Summary

Start from integers \(x \ge 2\) and \(y\), and define

$a_0=y,\qquad z_0=x,\qquad a_{t+1}=a_t^2 \bmod z_t,\qquad z_{t+1}=z_t+1.$

The sequence is counted until the first time \(a_t \le 1\). Since both \(0\) and \(1\) stay fixed under squaring, they are terminal values. Let \(\ell(x,y)\) be the number of terms from \(a_0\) through that first terminal term. Then

$g(x)=\max_y \ell(x,y),\qquad f(n)=\max_{2 \le x \le n} g(x).$

The direct search over all starting values is far too expensive, so the implementations reorganize the computation around the distinct residues that can still be alive after each step.

Mathematical Approach

The key observation is that the process quickly merges many different starts into the same residue. Once we track those shared residues instead of individual starts, the problem becomes manageable.

Step 1: Compress the Starting Values into First-Step Residues

For fixed \(x\), the first update depends only on \(y \bmod x\), so it is enough to consider residues \(0 \le y < x\). Moreover,

$ (x-y)^2 \equiv y^2 \pmod{x}, $

so the starts \(y\) and \(x-y\) merge immediately. For \(x>2\), the only starts worth keeping are \(2 \le y \le \lfloor x/2 \rfloor\), and only if the first residue is greater than \(1\). Define the first active frontier by

$A_1(x)=\left\{y^2 \bmod x : 2 \le y \le \left\lfloor \frac{x}{2} \right\rfloor,\ y^2 \bmod x > 1\right\}.$

If \(A_1(x)=\varnothing\), then every nontrivial start reaches \(0\) or \(1\) after one update, so \(g(x)=2\). The only special case is \(x=2\), where even the initial value is already forced to be terminal and \(g(2)=1\).

Step 2: Propagate the Entire Frontier at Once

For \(t \ge 1\), let \(A_t(x)\) be the set of residues still alive after \(t\) updates. The recurrence is

$A_{t+1}(x)=\left\{a^2 \bmod (x+t) : a \in A_t(x),\ a^2 \bmod (x+t) > 1\right\}.$

Each set update simultaneously represents every starting value whose trajectory has merged into one of those residues. When \(A_t(x)\) becomes empty, every trajectory has terminated, so the total length is \(t+1\). This replaces many separate simulations by one deduplicated state update per time step.

Step 3: Exploit Singleton Tails

If the frontier shrinks to a single residue, say

$A_t(x)=\{v\},$

then all surviving starting values now share exactly the same future. From that moment onward, there is no benefit in keeping a set representation: the rest is just the single trajectory starting from value \(v\) with current modulus \(x+t\). In other words, after coalescence the problem stops being a many-start search and becomes one ordinary modular-squaring chain.

Step 4: Derive an Interval Upper Bound for the Outer Maximum

To compute \(f(n)\), the implementations do not scan every \(x\) from \(2\) to \(n\). Instead they use a rigorous upper bound on intervals. Fix \(L \le u \le R\). Take any start \(y\) for the modulus \(u\), and advance the sequence for exactly \(R-u\) updates. At that point the current modulus is \(R\), and the current value is some residue \(b<R\). Therefore the remaining lifetime is at most \(g(R)\), which gives

$\ell(u,y) \le (R-u) + g(R).$

Taking the maximum over all \(y\) yields

$g(u) \le (R-u) + g(R) \le (R-L) + g(R).$

So every interval \([L,R]\) has the valid bound

$UB([L,R]) = g(R) + (R-L).$

If this bound is no better than the best value already known, the whole interval can be discarded safely.

Worked Example: \(x=10\)

The first distinct nonterminal residues are

$A_1(10)=\{2^2,3^2,4^2,5^2\}\bmod 10 \setminus \{0,1\}=\{4,5,6,9\}.$

Now propagate the frontier:

$A_2(10)=\{a^2 \bmod 11 : a \in A_1(10)\}\setminus\{0,1\}=\{3,4,5\},$

$A_3(10)=\{a^2 \bmod 12 : a \in A_2(10)\}\setminus\{0,1\}=\{4,9\},$

$A_4(10)=\{a^2 \bmod 13 : a \in A_3(10)\}\setminus\{0,1\}=\{3\}.$

At this point the frontier is a singleton, so only one tail remains. Continuing from value \(3\) with moduli \(14,15,16,\dots\) gives

$3 \to 9 \to 6 \to 4 \to 16 \to 4 \to 16 \to 16 \to \cdots \to 8 \to 0.$

The terminal \(0\) appears right after the step modulo \(32\), so the sequence has \(24\) terms in total and therefore

$g(10)=24.$

This example shows all the main ideas: first-step symmetry, deduplication of equal residues, and the switch to a single tail once every surviving start has merged.

How the Code Works

The C++, Python, and Java implementations follow the same logic. First they build the initial frontier by scanning only half of the possible starting residues, using the symmetry \(y \leftrightarrow x-y\), and discarding first-step residues \(0\) and \(1\). To avoid recomputing each square from scratch during that initial scan, the compiled implementations update consecutive squares incrementally.

While the frontier contains more than one residue, the implementation advances the whole set one modulus at a time. Large frontiers use a dense seen-array for fast deduplication; small frontiers use a sparse hash-based container to avoid touching every residue class of the modulus. Once the frontier becomes empty, the current length is the answer for that \(x\). Once the frontier becomes a singleton, the remaining steps are computed by direct modular squaring of that one value.

For the outer function \(f(n)\), the implementation first samples a coarse grid of \(x\)-values, caches the computed \(g(x)\), and places the intervals between sampled points into a max-priority queue keyed by the bound \(UB([L,R])\). It then repeatedly bisects only those intervals whose upper bound can still beat the current best answer. The C++ and Java implementations evaluate the initial grid points in parallel, while the Python implementation uses the same search strategy sequentially.

Complexity Analysis

For a fixed \(x\), building the first frontier costs \(O(x)\) time, because the scan goes through \(2,3,\dots,\lfloor x/2 \rfloor\). After that, one sparse frontier update costs expected \(O(|A_t(x)|)\), while one dense update costs \(O(x+t)\) because the deduplication array spans the whole current modulus. Therefore the running time for \(g(x)\) is output-sensitive: it is governed by the cumulative frontier sizes until extinction or singleton collapse, rather than by the number of possible starts alone.

The memory usage for one evaluation of \(g(x)\) is \(O(x+t_{\max})\) in dense phases and \(O(|A_t(x)|)\) in sparse phases. For \(f(n)\), if the interval search evaluates \(m\) different values of \(x\), the total cost is

$O\left(m + \sum_{i=1}^{m} T(x_i)\right),$

where \(T(x_i)\) is the cost of computing \(g(x_i)\). In the worst case \(m\) can still be \(O(n)\), but the interval bound usually prunes most subintervals long before that.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=693
2. Modular arithmetic: Wikipedia — Modular arithmetic
3. Quadratic residues: Wikipedia — Quadratic residue
4. Branch and bound: Wikipedia — Branch and bound
5. Priority queues: Wikipedia — Priority queue