Detailed mathematical approach

Problem Summary

Start from the reduced fraction \(1/n\). At each step replace \(x/y\) by \((x+1)/(y-1)\) and reduce again. If the first time the denominator becomes 1 the numerator is \(f(n)\), the problem asks for

$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$

A direct simulation for every cube is far too slow, so the key is to rewrite \(f(n)\) using the prime factors of \(n+1\).

Mathematical Approach

Write any reduced state as \(a/b\) with \(a,b>0\), and define its state sum by

$s=a+b.$

Initially we have \(a=1\), \(b=n\), hence \(s=n+1\).

One Step in Terms of the State Sum

From the reduced fraction \(a/b\), the unreduced next fraction is \((a+1)/(b-1)\). If

$d=\gcd(a+1,b-1),$

then after reduction

$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$

Since \(b=s-a\), we get

$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$

Therefore the next state sum is

$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$

So every reduction replaces the current sum \(s\) by a divisor of \(s\). Also, because \(a/b\) is reduced, \(\gcd(a,b)=1\), and thus \(\gcd(a,s)=1\).

A Phase Always Removes the Smallest Prime Factor

Suppose a phase starts at the special state \(1/(s-1)\). Let \(p\) be the smallest prime factor of \(s\). For every integer \(t\) with \(2 \le t \lt p\), we have \(\gcd(t,s)=1\). Hence no reduction happens while the numerator runs through \(1,2,\dots,p-1\).

When the unreduced numerator becomes \(p\), the reduction factor is

$d=\gcd(p,s)=p.$

The fraction at that moment is

$\frac{p}{s-p}\to \frac{1}{s/p-1}.$

So one full phase transforms the state sum by

$s\longmapsto \frac{s}{p},$

which means: a phase strips off the smallest prime factor of the current sum and resets the numerator to 1.

Closed Form for \(f(n)\)

Starting from \(s_0=n+1\), repeated phases remove the prime factors of \(n+1\) from smallest to largest. After all but the largest prime factor have been removed, the state sum becomes

$P^+(n+1)=\operatorname{LPF}(n+1),$

where \(\operatorname{LPF}\) denotes the largest prime factor.

If the current sum \(s\) is prime, then for every \(1\le a\le s-2\) we have \(\gcd(a+1,s)=1\), so no further reductions occur before the denominator reaches 1. Beginning from \(1/(s-1)\), the chain ends at \((s-1)/1\). Therefore

$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$

Worked Example: \(n=20\)

The iteration is

$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$

So \(f(20)=6\). Since \(20+1=21=3\cdot 7\), the closed form gives

$\operatorname{LPF}(21)-1=7-1=6,$

which matches the simulation exactly.

Cubic Specialization

For the required sum, set \(n=k^3\):

$f(k^3)=\operatorname{LPF}(k^3+1)-1.$

Using the factorization of a sum of cubes,

$k^3+1=(k+1)(k^2-k+1).$

If we define

$Q_k=k^2-k+1,$

then

$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$

This remains true even when the two factors share a factor 3, because the largest prime factor of a product is still the maximum of the largest prime factors of its factors. Hence

$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$

Sieving \(\operatorname{LPF}(k+1)\)

A standard largest-prime-factor sieve on \([2,L+1]\) fills an array \(\operatorname{LPF}(m)\) for every \(m\le L+1\). This immediately gives all values \(\operatorname{LPF}(k+1)\).

Sieving \(\operatorname{LPF}(Q_k)\)

For each \(k\), initialize

$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$

Now fix a prime \(p\). We need to know for which \(k\) the congruence

$Q_k\equiv0\pmod p$

holds. Multiplying by 4 gives

$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$

Therefore, for \(p>3\), roots exist exactly when \(-3\) is a quadratic residue modulo \(p\). If \(s^2\equiv -3\pmod p\), then

$k\equiv \frac{1\pm s}{2}\pmod p.$

The code finds the square root \(s\) with the Tonelli-Shanks algorithm. The prime \(p=3\) is special: the only class is \(k\equiv2\pmod 3\). The prime \(2\) never divides \(Q_k\), because \(k^2-k=k(k-1)\) is always even, so \(Q_k\) is always odd.

For every matching residue class, the code visits all corresponding \(k\), divides out all powers of \(p\) from \(\mathrm{rem}[k]\), and stores \(p\) as the current largest prime factor. After all primes \(p\le L\) are processed, any leftover \(\mathrm{rem}[k]>1\) must itself be prime and greater than \(L\); otherwise a composite remainder would have a prime divisor \(\le \sqrt{Q_k} \lt L\), which would already have been removed. That leftover is therefore the final large prime factor of \(Q_k\).

Worked Example: \(k=4\)

Here

$k^3+1=65=5\cdot 13,$

so

$f(64)=\operatorname{LPF}(65)-1=13-1=12.$

The decomposition gives \(k+1=5\) and \(Q_4=13\), so the larger prime factor is indeed \(13\).

How the Code Works

The implementation has two sieve layers. First it fills lpf_small for all integers up to \(L+1\). Then it stores every \(Q_k\) in the residual array rem and processes primes one by one. For each prime it computes the admissible residue classes of \(k\), strips that prime from all matching residuals, and updates lpf_q[k]. Finally it combines lpf_small[k+1] with lpf_q[k], subtracts 1, and accumulates the sum in a 128-bit total. The C++ version checks \(f(20)=6\), \(\sum_{k\le100} f(k^3)=118937\), and also matches a brute-force computation for small limits.

Complexity Analysis

The sieve for \(\operatorname{LPF}(k+1)\) is near-linear in practice. The \(Q_k\)-sieve visits about \(L/p\) or \(2L/p\) positions per relevant prime, so its dominant work is also near-linear up to the usual \(\sum_{p\le L}1/p\) factor, i.e. \(O(L\log\log L)\) with extra polylogarithmic cost for modular square roots. Memory usage is \(O(L)\) for the factor tables and residual buffer.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=343
2. Largest prime factor: Wikipedia - Prime factor
3. Tonelli-Shanks algorithm: Wikipedia - Tonelli-Shanks algorithm
4. Modular square roots: Wikipedia - Quadratic residue
5. Sum of cubes factorization: Wikipedia - Polynomial factorization identities