Detailed mathematical approach

Problem Summary

Define

$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$

The target quantity is the gcd-weighted double sum

$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$

with the Project Euler input \(N=10^8\). A direct evaluation is impossible: the outer expression ranges over \(N^2\) ordered pairs, while each weight \(w_n\) is itself a nested lookup inside the same cubic recurrence. The solution works only because the gcd part can be counted arithmetically and the recurrence part becomes periodic.

Mathematical Approach

The derivation has two main pieces. First, the gcd-weighted double sum is reorganized by possible gcd values and reduced to the summatory totient function. Second, the cubic orbit is shown to have a short preperiod and a manageable cycle, which turns the nested weight \(w_n=x_{x_n}\) into a table lookup inside one periodic tail.

Grouping Ordered Pairs by Their GCD

For a fixed ordered pair \((i,j)\), write

$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$

Once \(d\) is fixed, the reduced pair \((a,b)\) must satisfy

$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor.$

If we define

$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\},$

then every ordered pair whose gcd is exactly \(d\) contributes the same weight \(w_d\). Therefore

$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$

This is the first decisive simplification: the original two-dimensional sum is replaced by a one-dimensional sum over gcd classes.

Counting Coprime Pairs with the Summatory Totient

The quantity \(R(q)\) counts ordered coprime pairs inside a \(q\times q\) grid. For each denominator \(m\ge 2\), there are exactly \(\varphi(m)\) reduced numerators \(1\le a<m\). Because the pairs are ordered, both \((a,m)\) and \((m,a)\) appear, and the diagonal contributes the extra point \((1,1)\). Hence

$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$

where

$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$

Substituting into the gcd decomposition gives the exact formula used by the implementations:

$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$

A Recurrence for \(\Phi(n)\)

The remaining arithmetic bottleneck is the summatory totient \(\Phi(n)\). The standard identity behind the code is

$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$

Isolating the \(k=1\) term yields

$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$

The quotient \(\left\lfloor n/k\right\rfloor\) stays constant over intervals. If \(q=\left\lfloor n/\ell\right\rfloor\), then the last index with the same quotient is

$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$

Grouping equal quotients gives the block form

$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$

This is what makes large values of \(\Phi\) practical: the recursion touches only the distinct quotient values, not all integers one by one.

The Periodic Tail of the Cubic Orbit

The recurrence

$f(t)\equiv (t-1)^3+2 \pmod{M}$

acts on a finite state space, so the orbit starting at \(x_1=1\) must eventually repeat. Floyd's cycle-finding method, applied to the actual orbit, gives a preperiod of length \(53\) and a cycle of length \(33705\). Therefore the periodic tail begins at \(x_{54}\), and for every \(r\ge 54\),

$x_r=x_{\,54+((r-54)\bmod 33705)}.$

The implementations therefore precompute the orbit through \(x_{33759}\), which is enough to cover the whole non-periodic prefix and every periodic representative that is later queried.

Turning \(w_n=x_{x_n}\) into a Tail Lookup

The nested index is the subtle part of the problem. Since the same cubic polynomial defines the sequence, reducing the state modulo the cycle length \(33705\) is compatible with iteration. Define the companion orbit

$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$

Then \(y_n\equiv x_n\pmod{33705}\) for every \(n\). On the actual orbit, the only inner indices that need the non-periodic prefix are the first four:

$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$

After these exceptional cases, the inner index is lifted back into the periodic tail by

$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}.$

So the huge-looking object \(x_{x_n}\) is never evaluated by expanding the sequence out to position \(x_n\); it is recovered from one residue modulo the cycle length and one lookup inside the stored tail.

Compressing the Outer Sum by Equal Quotients

The outer factor depends on \(d\) only through

$q=\left\lfloor \frac{N}{d}\right\rfloor.$

If this quotient is constant on an interval \(d\in [L,R]\), then the whole block contributes

$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$

This is the second major compression. Instead of multiplying each \(w_d\) by its totient factor separately, the implementations maintain a running sum of the weights over each quotient block and apply the factor once. Since \(\left\lfloor N/d\right\rfloor\) has only \(O(\sqrt{N})\) distinct values, the expensive arithmetic attached to \(\Phi\) is heavily reused.

Worked Example: \(N=4\)

The first terms of the orbit are

$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$

Hence the first gcd-weights are

$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$

For \(N=4\), the relevant summatory-totient values are

$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$

so

$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$

The divisor formula gives

$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$

Numerically,

$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$

This tiny example already shows the block structure: \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\), so the last two divisor classes share the same factor.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical plan. They never iterate over all \(N^2\) pairs. Instead, they combine orbit preprocessing, a fast summatory-totient routine, and a quotient-block sweep over the divisor index \(d\).

Orbit Preprocessing

The implementation first applies Floyd's tortoise-and-hare method to the cubic recurrence and recovers the exact preperiod length \(53\) and cycle length \(33705\). It then stores the actual orbit values through the end of the useful periodic window, so every later access to \(x_k\) inside the tail becomes an array lookup. In parallel with the main divisor sweep, it also advances the same recurrence modulo \(33705\) so that the residue class of the inner index is always available.

Totient Preprocessing and Memoized \(\Phi\)

For smaller arguments, the implementation builds a prefix table of \(\varphi(n)\) up to \(10^6\) with a linear sieve. For larger arguments, it uses the quotient-block recurrence for \(\Phi(n)\) together with memoization. This hybrid is important: values below the cutoff are answered in constant time from the prefix table, while larger values reuse the same cached quotient results again and again.

Blockwise Accumulation of the Answer

The divisor range \(1\le d\le N\) is partitioned into maximal intervals on which \(\left\lfloor N/d\right\rfloor\) is constant. While sweeping those intervals from left to right, the implementation generates the next weights \(w_d\), updates a running prefix sum of weights modulo \(M\), extracts the segment sum for the current block, computes the factor \(2\Phi(q)-1\), and adds the block contribution. This is why the final loop is linear in \(N\) rather than quadratic in \(N\).

Complexity Analysis

Let \(L=\min(N,10^6)\). Cycle detection and orbit precomputation take \(O(53+33705)\) time and \(O(33705)\) memory, which are negligible compared with the main sweep. Building the totient prefix table costs \(O(L)\) time and \(O(L)\) memory because the sieve is linear. The memoized summatory-totient recurrence visits only \(O(\sqrt{N})\) distinct quotient values, and the quotient partition itself also has \(O(\sqrt{N})\) blocks.

However, the algorithm still generates every weight \(w_d\) for \(1\le d\le N\), so the dominant running time is \(O(N)\). Overall, the implementations run in \(O(N+L+\sqrt{N})\) time and use \(O(L+\sqrt{N}+33705)\) memory.

Footnotes and References

1. Project Euler Problem 915: https://projecteuler.net/problem=915
2. Euler's totient function: Wikipedia - Euler's totient function
3. Totient summatory function: Wikipedia - Totient summatory function
4. Cycle detection: Wikipedia - Cycle detection
5. Coprime integers and visible lattice points: Wikipedia - Coprime integers