Problem 343: Fractional Sequences

Project Euler Problem 343 Solution

EulerSolve provides an optimized solution for Project Euler Problem 343, Fractional Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

Detailed mathematical approach

Problem Summary

Fix a positive integer $k$. The first fraction of a finite list is $a_1=1/k$, and every later term comes from its predecessor $x_{i-1}/y_{i-1}$ as

$$a_i=\frac{x_{i-1}+1}{y_{i-1}-1},$$

reduced to lowest terms: add one to the numerator, subtract one from the denominator, cancel. The list stops once the reduced denominator is $y_i=1$, so the fraction has become an integer $n$, and we set $f(k)=n$. The statement works out $k=20$ in full,

$$\begin{aligned}\tfrac1{20}\to\tfrac2{19}\to\tfrac3{18}&=\tfrac16\to\tfrac25\\&\to\tfrac34\to\tfrac43\to\tfrac52\to\tfrac61,\end{aligned}$$

so $f(20)=6$. It also gives $f(1)=1$, $f(2)=2$ and $f(3)=1$, with the checkpoint $\sum f(k^3)=118937$ for $1\le k\le100$. We must evaluate ^[1]

$$\sum_{k=1}^{L}f\bigl(k^3\bigr),\qquad L=2\cdot10^{6}.$$

The whole solution in one sentence: $f(k)$ is exactly one less than the largest prime factor of $k+1$, so the task becomes a sieve that finds the largest prime factor of $(k+1)(k^2-k+1)$ for every $k$ up to two million.

Mathematical Approach

Simulating the list is out of the question: one term of the sum alone can need more than $10^{12}$ steps. Everything follows instead from a single observation — the numerator and the denominator add up to a number that can only shrink, and only by whole prime factors.

1. The sum of numerator and denominator

Let $a/b$ be one of the reduced fractions and call $s=a+b$ its state sum. At the start $a=1$ and $b=k$, so $s=k+1$. While $b\gt1$ the list continues, and the numerator and denominator of the unreduced successor $(a+1)/(b-1)$ again add up to $s$, so their common factor is

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

because $b-1=s-(a+1)$ and $\gcd(u,v)=\gcd(u,u+v)$. Cancelling gives $a'=(a+1)/d$ and $b'=(b-1)/d$, so the new state sum is

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

So each state sum divides the previous one, and every state sum divides $k+1$. Also, $a/b$ being reduced means $\gcd(a,b)=1$, which by the same identity is $\gcd(a,s)=1$: the numerator of a reduced state is always coprime to its state sum.

2. Nothing cancels until the numerator reaches the smallest prime factor

Call a state fresh when its numerator is $1$, so that it looks like $1/(s-1)$; the start is fresh. Let $p$ be the smallest prime factor of $s$, and suppose for the moment that $s$ is not prime. Each numerator $2,3,\ldots,p-1$ is a product of primes smaller than $p$, hence coprime to $s$, so $d=1$ and the step merely shifts one unit from the denominator to the numerator:

$$\frac{1}{s-1}\to\frac{2}{s-2}\to\cdots\to\frac{p-1}{s-p+1}.$$

The next step produces $p/(s-p)$, where $d=\gcd(p,s)=p$, and cancelling gives the fresh state $1/(s/p-1)$. One such phase therefore costs $p-1$ steps and replaces $s$ by $s/p$. The denominator stays positive, since $s$ composite forces $p\le\sqrt s$. For $k=64$, with $s=65$ and $p=5$, the chain runs $1/64\to2/63\to3/62\to4/61$ without cancelling, and the fourth step turns $5/60$ into $1/12$.

3. The closed form for $f(k)$

If $s$ is prime, nothing can cancel at all: every numerator produced along the way satisfies $2\le a+1\le s-1$, so $\gcd(a+1,s)=1$. The chain from $1/(s-1)$ then walks all the way to $(s-1)/1$ and stops there with value $s-1$.

Starting from $s=k+1$, Section 2 strips the smallest prime factor again and again and stops as soon as a single prime is left. By unique factorisation that survivor is the largest prime factor $P(k+1)$. ^[2] Hence

$$f(k)=P(k+1)-1,$$

with the usual convention $P(1)=1$. ^[3] The statement's three values come out right: $P(2)-1=1$, $P(3)-1=2$ and $P(4)-1=1$. For $k=1$ the sum $2$ is already prime and $1/1$ is already an integer, so the list stops at once.

4. Worked example: the chain for $k=29$

Here $s=30=2\cdot3\cdot5$, so the state sums must run $30\to15\to5$. The first phase has $p=2$ and cancels immediately; the second has $p=3$ and takes one quiet step first; the third sum is prime and runs to the end:

$$\begin{aligned}\tfrac1{29}\to\tfrac1{14}&\to\tfrac2{13}\to\tfrac14\\&\to\tfrac23\to\tfrac32\to\tfrac41,\end{aligned}$$

so $f(29)=4=P(30)-1$. Equal prime factors are stripped one at a time: $k=124$ has $s=125=5^3$ with sums $125\to25\to5$, so $f(124)=4$ as well.

Counting phases also gives the length of the list. A phase for $p$ costs $p-1$ steps and the closing prime phase costs $P-2$, so if $k+1=\prod p_i^{e_i}$ the number of steps is $\sum e_i(p_i-1)-1$. For $k=20$ that is $2+6-1=7$ steps, hence eight fractions, exactly as printed; for $k=29$ it is $1+2+4-1=6$ steps and seven fractions.

5. Why a direct simulation is hopeless

The sum needs $f(k^3)$, so the state sum starts at $k^3+1$, about $8\cdot10^{18}$ at the top of the range. Take $k=1999999$: here $k+1=2^7\cdot5^6$ and $k^2-k+1=3999994000003$ is prime, so the step count above is $3999994000032$ — roughly $4\cdot10^{12}$ additions and reductions for one of two million terms.

6. Cubes split into two factors

Multiplying out $(k+1)(k^2-k+1)$${}=k^3-k^2+k+k^2-k+1$ gives

$$k^3+1=(k+1)\bigl(k^2-k+1\bigr),$$

and we write $Q_k=k^2-k+1$. ^[4] Every prime factor of a product divides one of the factors, so

$$P\bigl(k^3+1\bigr)=\max\bigl(P(k+1),P(Q_k)\bigr)$$

and $f(k^3)=\max(P(k+1),P(Q_k))-1$. The two factors are almost coprime: $Q_k=(k+1)(k-2)+3$, so $\gcd(k+1,Q_k)=\gcd(k+1,3)$, which is $1$ or $3$. An overlap does no harm, since a maximum of largest prime factors is correct whether or not the factors share something. For $k=1$ the second factor is $Q_1=1$, and the convention $P(1)=1$ gives $f(1)=\max(2,1)-1=1$.

7. Which primes divide $Q_k$

The factor $k+1$ is at most $L+1$, and an ordinary sieve handles it. The factor $Q_k$ reaches $3999998000001$, far too large to factorise two million times — but very few primes divide it at all. Below we work modulo $p$: $u\equiv v$ means $p$ divides $u-v$, and the multiplicative order of $u$ is the least $d\gt0$ with $u^d\equiv1$, which divides every exponent that yields $1$. ^[5]

Never $p=2$. $k^2-k=k(k-1)$ is a product of consecutive integers, hence even, so $Q_k$ is odd.
$p=3$ exactly when $k\equiv2\pmod 3$. Writing $k=3m+2$ gives $Q_k=3(3m^2+3m+1)$, and the bracket leaves remainder $1$ modulo $3$, so $3$ divides $Q_k$ once and $9$ never does.
$p\gt3$ only when $p\equiv1\pmod 6$. If $p\mid Q_k$ then $p\mid k^3+1$, so $k^3\equiv-1$ and $k^6\equiv1$, and the order of $k$ divides $6$. Order $1$ or $3$ would give $1\equiv-1$ and hence $p=2$; order $2$ means $k\equiv-1$, which makes $Q_k\equiv3$; so the order is $6$, and since it divides $p-1$ we get $6\mid p-1$.

Of the $148933$ primes up to $L$, only $74412$ satisfy $p\equiv1\pmod 3$, so together with $p=3$ just $74413$ primes can do any work.

8. Finding the k that a prime hits

For a fixed prime $p\gt3$ we want every $k$ with $p\mid Q_k$. Multiplying by $4$ completes the square,

$$4Q_k=(2k-1)^2+3,$$

so $p\mid Q_k$ is the same as $(2k-1)^2\equiv-3\pmod p$. Such a root exists only if $-3$ is a quadratic residue, which the code tests with Euler's criterion, $(-3)^{(p-1)/2}\equiv1\pmod p$. ^[6] That test makes skipping a prime safe: if $-3$ is a non-residue, no $k$ at all has $p\mid Q_k$. By Section 7 a root needs $p\equiv1\pmod 3$, and at $L=2\cdot10^{6}$ exactly those $74412$ primes pass.

Given one square root $t$, the two residue classes are

$$k\equiv\frac{1\pm t}{2}\pmod p,$$

where $1/2$ stands for the inverse $(p+1)/2$. The root itself comes from the Tonelli–Shanks algorithm. ^[7] For $p\equiv3\pmod 4$ one exponentiation suffices, since $(-3)^{(p+1)/2}=(-3)\cdot(-3)^{(p-1)/2}\equiv-3$. For $p=7$, say, $-3\equiv4$ and $t=2$, giving $k\equiv3$ and $k\equiv5$ modulo $7$ — and indeed $Q_3=7$, $Q_5=21=3\cdot7$ and $Q_{10}=91=7\cdot13$.

9. Two sieves and the large leftover

The first sieve walks the primes $p\le L+1$ in increasing order and writes $p$ into every multiple of $p$. The last write to a position wins, so each entry ends up holding its own largest prime factor, which gives $P(k+1)$ for every $k$.

The second sieve starts with $Q_k$ in a residual array. For each of the $74413$ relevant primes it computes the classes of Section 8 and, at every position in them, divides all powers of $p$ out of the residual and records $p$. As the primes come in increasing order, the last one recorded is the largest prime factor of $Q_k$ up to $L$ — or the initial $1$, if there is none.

A residual still greater than $1$ must then be a single prime larger than $L$: all its prime factors exceed $L$, and a product of two would be at least $(L+1)^2$, whereas $Q_k\le L^2-L+1\lt L^2$. At $L=2\cdot10^{6}$ this happens for $1373993$ values of $k$, almost $69\%$, so the large leftover is the normal outcome.

10. Worked example: the sum up to $k=10$

With both tables in hand every term is a maximum of two numbers. For $k=1,\ldots,10$ the factor $k+1$ supplies the largest prime factors $2,3,2,5,3,7,2,3,5,11$, while $Q_k=1,3,7,13,21,31,43,57,73,91$ supplies $1,3,7,13,7,31,43,19,73,13$. Taking maxima and subtracting one gives

$$1,2,6,12,6,30,42,18,72,12,$$

whose sum is $201$. The quadratic factor usually decides the value, but not always: at $k=82$ it is $Q_{82}=6643=7\cdot13\cdot73$ while $k+1=83$ is prime, so $f(82^3)=82$.

How the Code Works

The C++ program ^[9] keeps its parameters in Options, which parse_arguments fills from --limit= and --skip-checkpoints. solve(limit) then carries out Section 9. It first builds lpf_small, the largest prime factors on $[2,\,\mathtt{limit}+1]$, collecting the primes up to limit in primes as it goes. Then rem[k] is filled with $Q_k$ and lpf_q[k] with $1$, the convention $P(1)=1$ of Section 6.

The prime loop skips $p=2$ (Section 7) and handles $p=3$ by its single class roots = {2}. For $p\gt3$ it applies Euler's criterion through pow_mod, skips the prime when $-3$ is a non-residue, calls tonelli_shanks and forms the two classes with inv2. Modular products go through mul_mod, which widens to __uint128_t. The inner loop then divides out every power of $p$ and stores $p$ in lpf_q; its test value % p != 0 is defensive only, since a position in the class satisfies $p\mid Q_k$ and only other primes have been removed. A final loop lifts any leftover residual into lpf_q, takes the maximum with lpf_small[k + 1], subtracts one and accumulates in a 128-bit total printed by to_string_u128.

The checks live in run_checkpoints. brute_f_of_k is the literal definition — increment the numerator, decrement the denominator, divide by std::gcd — and it confirms $f(20)=6$. Then solve(100) must return $118937$, and solve(300) must equal brute_sum_k_cubed(300), which simulates all $300$ chains. That helper computes k * k * k in int, so it is usable only while $k^3$ fits in 32 bits, up to $k=1290$; the checkpoint calls it with $300$.

The ports carry the same algorithm. Python's solver is solve_343(limit), wrapped by a zero-argument solve() that fixes the limit; mul_mod, pow_mod and tonelli_shanks mirror their C++ counterparts, with pow_mod delegating to the built-in three-argument pow, and arbitrary-precision integers remove every width question. Java wraps solve(int limit) the same way, keeps int[] and long[] arrays, accumulates in a plain long as the bound below justifies, and routes mulMod through BigInteger although the products stay below $4\cdot10^{12}$. Neither port has the argument parser or the checkpoints.

Complexity and Verification

The first sieve writes $\sum_{p\le L+1}\lfloor (L+1)/p\rfloor$ values, which is $5808802$ at the target, about $2.90$ per position. Mertens' second theorem gives $\sum_{p\le x}1/p=\ln\ln x+M+o(1)$, so this is $\Theta(L\log\log L)$. ^[8]

The second sieve sweeps $148825$ residue classes belonging to $74413$ primes and visits $4588668$ positions, about $2.29$ per value of $k$. Each prime contributes at most two classes of at most $L/p+1$ positions, so the same theorem bounds the visits by $O(L\log\log L)$. A visit costs one division plus the rare extra ones for higher powers of $p$, at most $2\log_2 L$ per position because $Q_k\lt L^2$.

Root finding is counted separately. Every prime $p\gt3$ pays one Euler-criterion exponentiation, $148931$ in all, and the $74412$ that pass call tonelli_shanks, which repeats that test; $37232$ of them have $p\equiv3\pmod 4$ and take the shortcut, while $37180$ enter the general loop. That loop first hunts for a quadratic non-residue by trying $2,3,4,\ldots$, a search with no useful worst-case bound, though here it tested $140033$ candidates in all and never needed a value above $41$; the rest of a call costs $O(\log^2 p)$. Neither part dominates at this size: at roughly $2\log_2 p$ modular multiplications each, the Euler tests alone come to some $6\cdot10^{6}$ against $4.6\cdot10^{6}$ sieve visits.

Memory is $O(L)$: one 32-bit array of $8.0$ MB and two 64-bit arrays of $16.0$ MB each, $40$ MB in all, while Python holds references to integer objects and needs more. Measured with GCC 16 at -O2, OpenJDK 17 and CPython 3.11 on one machine, peak resident memory was $44.4$ MB for C++ and $216.3$ MB for Python, and the running times were $0.11$ s for C++, $0.19$ s with the checkpoints, $0.53$ s for Java and $1.8$ s for Python.

Nothing overflows. Each $Q_k$ is at most $3999998000001$, well inside 64 bits. For $k\ge2$ we have $P(k^3+1)\le\max(k+1,Q_k)=Q_k$, so $f(k^3)\le k^2-k$, and the whole sum is at most

$$1+\sum_{k=2}^{L}\bigl(k^2-k\bigr)=2666666666666000001,$$

under $29\%$ of $2^{63}-1=9223372036854775807$. Java's signed long is therefore safe, and the C++ 128-bit accumulator is simply conservative.

For this article the method was checked independently of the packaged checkpoints. Direct simulation confirms $f(k)=P(k+1)-1$ for every $k\le3000$, and the fast method agrees with a full simulation of all cubes at $L=1,2,3,5,10,50,100,200,300$, giving $1$, $3$, $9$, $27$, $201$, $11835$, $118937$, $866597$ and $2333647$. Factorising $k^3+1$ independently for every $k\le50000$ returns $5670713547181$, matching the sieve. At the full limit the C++, Python and Java programs produced identical output, as did a separate reimplementation whose residue classes come from a different modular-square-root routine. The structural claims hold too: for $k\le200000$, $3$ divides $Q_k$ exactly once when $k\equiv2\pmod 3$ and never otherwise; for limits up to $20000$ every leftover residual is prime and larger than $L$; and the defensive division test never fires.

Footnotes and References

The references below supply the background facts used above. The state-sum argument, the closed form for $f(k)$, the step count and the operation counts are derived explicitly in this article.

Project Euler 343 — Fractional Sequences. The official statement defines a₁ = 1/k and aᵢ = (xᵢ₋₁ + 1)/(yᵢ₋₁ − 1) reduced to lowest terms, stops when the denominator is 1 and sets f(k) to that integer. It works out k = 20 ending at 6/1, gives f(1) = 1, f(2) = 2, f(3) = 1 and the sum of f(k³) as 118937 for 1 ≤ k ≤ 100, then asks for that sum over 1 ≤ k ≤ 2 × 10⁶.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes, unique apart from their order, and 1 is covered by the empty-product convention. Uniqueness makes "the largest prime factor" well defined in Section 3 and the maximum rule of Section 6 valid.
OEIS A006530 — greatest prime dividing n. The greatest prime factor of n for n ≥ 2, with the purely conventional a(1) = 1; it begins 1, 2, 3, 2, 5, 3, 7, 2, 3, 5, 11. This is the function P, and the convention the code stores in lpf_q[1].
OEIS A002061 — central polygonal numbers. The sequence a(n) = n² − n + 1 is the factor Q of Section 6, beginning 1, 1, 3, 7, 13, 21, 31, 43, 57, 73, 91. The entry notes that all terms are odd and that 3 divides a(3k − 1), and identifies the prime factors with OEIS A007645: 3 together with the primes 6m + 1, the classification proved in Section 7.
Wikipedia — Multiplicative order. For a coprime to n the order is the smallest positive k with aᵏ ≡ 1 (mod n), and by Lagrange's theorem it divides φ(n). Since φ(p) = p − 1, that is what turns order 6 into 6 | p − 1.
Wikipedia — Euler's criterion. For an odd prime p and a coprime to p, a raised to (p − 1)/2 is congruent to 1 modulo p when a is a square modulo p and to −1 otherwise. This is the test the programs run before seeking a square root of −3.
Wikipedia — Tonelli–Shanks algorithm. It solves r² ≡ n modulo a prime p by writing p − 1 = Q·2ˢ with Q odd, searching for a non-residue z and iterating until the error factor is 1. It also gives the shortcut used here for p ≡ 3 (mod 4), with candidates ±n^((p+1)/4).
Wikipedia — Mertens' theorems. The second theorem states that the sum of 1/p over the primes p up to n equals log log n plus the Meissel–Mertens constant M plus a term tending to 0. It supplies the log log factor in both sieve bounds; the article names M without giving its value.
C++ — Euler343.cpp. The linked immutable revision contains solve, tonelli_shanks, brute_f_of_k, brute_sum_k_cubed and run_checkpoints. It is the implementation record for this article; the extra comparisons above are not packaged in it.

Problem 343 source code

C++

#include <algorithm>
#include <cstdint>
#include <functional>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Options {
    int limit = 2'000'000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    u64 result = 1ULL;
    base %= mod;
    while (exp > 0ULL) {
        if (exp & 1ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

u64 tonelli_shanks(const u64 n, const u64 p) {
    if (n == 0ULL) {
        return 0ULL;
    }
    if (p == 2ULL) {
        return n;
    }
    if (pow_mod(n, (p - 1ULL) / 2ULL, p) != 1ULL) {
        return 0ULL;
    }
    if ((p & 3ULL) == 3ULL) {
        return pow_mod(n, (p + 1ULL) / 4ULL, p);
    }

    u64 q = p - 1ULL;
    int s = 0;
    while ((q & 1ULL) == 0ULL) {
        q >>= 1ULL;
        ++s;
    }

    u64 z = 2ULL;
    while (pow_mod(z, (p - 1ULL) / 2ULL, p) != p - 1ULL) {
        ++z;
    }

    u64 c = pow_mod(z, q, p);
    u64 x = pow_mod(n, (q + 1ULL) / 2ULL, p);
    u64 t = pow_mod(n, q, p);
    int m = s;

    while (t != 1ULL) {
        int i = 1;
        u64 t2i = mul_mod(t, t, p);
        while (i < m && t2i != 1ULL) {
            t2i = mul_mod(t2i, t2i, p);
            ++i;
        }
        const u64 b = pow_mod(c, 1ULL << (m - i - 1), p);
        x = mul_mod(x, b, p);
        t = mul_mod(t, mul_mod(b, b, p), p);
        c = mul_mod(b, b, p);
        m = i;
    }
    return x;
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string s;
    while (value > 0) {
        const int d = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + d));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 solve(const int limit) {
    std::vector<int> lpf_small(static_cast<std::size_t>(limit + 2), 0);
    std::vector<int> primes;
    for (int i = 2; i <= limit + 1; ++i) {
        if (lpf_small[static_cast<std::size_t>(i)] == 0) {
            for (int j = i; j <= limit + 1; j += i) {
                lpf_small[static_cast<std::size_t>(j)] = i;
            }
            if (i <= limit) {
                primes.push_back(i);
            }
        }
    }

    std::vector<u64> rem(static_cast<std::size_t>(limit + 1), 0ULL);
    std::vector<u64> lpf_q(static_cast<std::size_t>(limit + 1), 1ULL);
    for (int k = 1; k <= limit; ++k) {
        const u64 kk = static_cast<u64>(k);
        rem[static_cast<std::size_t>(k)] = kk * kk - kk + 1ULL;
    }

    for (int p : primes) {
        if (p == 2) {
            continue;
        }
        std::vector<int> roots;
        if (p == 3) {
            roots.push_back(2);
        } else {
            const u64 mod = static_cast<u64>(p);
            const u64 neg_three = (mod + mod - 3ULL) % mod;
            if (pow_mod(neg_three, (mod - 1ULL) / 2ULL, mod) != 1ULL) {
                continue;
            }
            const u64 sq = tonelli_shanks(neg_three, mod);
            const u64 inv2 = (mod + 1ULL) / 2ULL;
            const int r1 = static_cast<int>(mul_mod((1ULL + sq) % mod, inv2, mod));
            const int r2 = static_cast<int>(mul_mod((1ULL + mod - sq) % mod, inv2, mod));
            roots.push_back(r1);
            if (r2 != r1) {
                roots.push_back(r2);
            }
        }

        for (int root : roots) {
            int k = root;
            if (k == 0) {
                k += p;
            }
            for (; k <= limit; k += p) {
                u64& value = rem[static_cast<std::size_t>(k)];
                if (value % static_cast<u64>(p) != 0ULL) {
                    continue;
                }
                while (value % static_cast<u64>(p) == 0ULL) {
                    value /= static_cast<u64>(p);
                }
                lpf_q[static_cast<std::size_t>(k)] = static_cast<u64>(p);
            }
        }
    }

    u128 total = 0;
    for (int k = 1; k <= limit; ++k) {
        const u64 rem_val = rem[static_cast<std::size_t>(k)];
        if (rem_val > 1ULL) {
            lpf_q[static_cast<std::size_t>(k)] = std::max(lpf_q[static_cast<std::size_t>(k)], rem_val);
        }
        const u64 lpf_kp1 = static_cast<u64>(lpf_small[static_cast<std::size_t>(k + 1)]);
        const u64 largest = std::max(lpf_kp1, lpf_q[static_cast<std::size_t>(k)]);
        total += static_cast<u128>(largest - 1ULL);
    }
    return total;
}

u64 brute_f_of_k(const int k) {
    u64 x = 1ULL;
    u64 y = static_cast<u64>(k);
    while (y != 1ULL) {
        ++x;
        --y;
        const u64 g = std::gcd(x, y);
        x /= g;
        y /= g;
    }
    return x;
}

u128 brute_sum_k_cubed(const int limit) {
    u128 s = 0;
    for (int k = 1; k <= limit; ++k) {
        const int kc = k * k * k;
        s += static_cast<u128>(brute_f_of_k(kc));
    }
    return s;
}

bool run_checkpoints() {
    if (brute_f_of_k(20) != 6ULL) {
        std::cerr << "Checkpoint failed: f(20)=6" << '\n';
        return false;
    }
    if (to_string_u128(solve(100)) != "118937") {
        std::cerr << "Checkpoint failed: sum_{k<=100} f(k^3)=118937" << '\n';
        return false;
    }
    if (solve(300) != brute_sum_k_cubed(300)) {
        std::cerr << "Checkpoint failed: brute-force cross-check for limit=300" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << to_string_u128(solve(options.limit)) << '\n';
    return 0;
}

Python

import math

def mul_mod(a, b, mod):
    return (a * b) % mod

def pow_mod(base, exp, mod):
    return pow(base, exp, mod)

def tonelli_shanks(n, p):
    if n == 0:
        return 0
    if p == 2:
        return n
    if pow_mod(n, (p - 1) // 2, p) != 1:
        return 0
    if (p % 4) == 3:
        return pow_mod(n, (p + 1) // 4, p)
        
    q = p - 1
    s = 0
    while (q % 2) == 0:
        q //= 2
        s += 1
        
    z = 2
    while pow_mod(z, (p - 1) // 2, p) != p - 1:
        z += 1
        
    c = pow_mod(z, q, p)
    x = pow_mod(n, (q + 1) // 2, p)
    t = pow_mod(n, q, p)
    m = s
    
    while t != 1:
        i = 1
        t2i = mul_mod(t, t, p)
        while i < m and t2i != 1:
            t2i = mul_mod(t2i, t2i, p)
            i += 1
            
        b = pow_mod(c, 1 << (m - i - 1), p)
        x = mul_mod(x, b, p)
        t = mul_mod(t, mul_mod(b, b, p), p)
        c = mul_mod(b, b, p)
        m = i
        
    return x

def solve_343(limit):
    lpf_small = [0] * (limit + 2)
    primes = []
    
    for i in range(2, limit + 2):
        if lpf_small[i] == 0:
            for j in range(i, limit + 2, i):
                lpf_small[j] = i
            if i <= limit:
                primes.append(i)
                
    rem = [0] * (limit + 1)
    lpf_q = [1] * (limit + 1)
    
    for k in range(1, limit + 1):
        rem[k] = k * k - k + 1
        
    for p in primes:
        if p == 2:
            continue
        roots = []
        if p == 3:
            roots.append(2)
        else:
            neg_three = (2 * p - 3) % p
            if pow_mod(neg_three, (p - 1) // 2, p) != 1:
                continue
            sq = tonelli_shanks(neg_three, p)
            inv2 = (p + 1) // 2
            r1 = ((1 + sq) % p * inv2) % p
            r2 = ((1 + p - sq) % p * inv2) % p
            roots.append(r1)
            if r2 != r1:
                roots.append(r2)
                
        for root in roots:
            k = root
            if k == 0:
                k += p
            for j in range(k, limit + 1, p):
                val = rem[j]
                if val % p != 0:
                    continue
                while val % p == 0:
                    val //= p
                rem[j] = val
                lpf_q[j] = p
                
    total = 0
    for k in range(1, limit + 1):
        rem_val = rem[k]
        if rem_val > 1:
            lpf_q[k] = max(lpf_q[k], rem_val)
        lpf_kp1 = lpf_small[k + 1]
        largest = max(lpf_kp1, lpf_q[k])
        total += largest - 1
        
    return total

def solve():
    ans = solve_343(2000000)
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler343 {

    static long mulMod(long a, long b, long mod) {
        return (long) ((((java.math.BigInteger.valueOf(a)).multiply(java.math.BigInteger.valueOf(b)))
                .remainder(java.math.BigInteger.valueOf(mod))).longValue());
    }

    static long powMod(long base, long exp, long mod) {
        long result = 1;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) == 1) {
                result = mulMod(result, base, mod);
            }
            base = mulMod(base, base, mod);
            exp >>= 1;
        }
        return result;
    }

    static long tonelliShanks(long n, long p) {
        if (n == 0)
            return 0;
        if (p == 2)
            return n;
        if (powMod(n, (p - 1) / 2, p) != 1)
            return 0;
        if ((p % 4) == 3)
            return powMod(n, (p + 1) / 4, p);

        long q = p - 1;
        int s = 0;
        while ((q % 2) == 0) {
            q /= 2;
            s++;
        }

        long z = 2;
        while (powMod(z, (p - 1) / 2, p) != p - 1) {
            z++;
        }

        long c = powMod(z, q, p);
        long x = powMod(n, (q + 1) / 2, p);
        long t = powMod(n, q, p);
        int m = s;

        while (t != 1) {
            int i = 1;
            long t2i = mulMod(t, t, p);
            while (i < m && t2i != 1) {
                t2i = mulMod(t2i, t2i, p);
                i++;
            }
            long b = powMod(c, 1L << (m - i - 1), p);
            x = mulMod(x, b, p);
            t = mulMod(t, mulMod(b, b, p), p);
            c = mulMod(b, b, p);
            m = i;
        }
        return x;
    }

    static long solve(int limit) {
        int[] lpfSmall = new int[limit + 2];
        List<Integer> primes = new ArrayList<>();

        for (int i = 2; i <= limit + 1; i++) {
            if (lpfSmall[i] == 0) {
                for (int j = i; j <= limit + 1; j += i) {
                    lpfSmall[j] = i;
                }
                if (i <= limit) {
                    primes.add(i);
                }
            }
        }

        long[] rem = new long[limit + 1];
        long[] lpfQ = new long[limit + 1];
        Arrays.fill(lpfQ, 1);

        for (int k = 1; k <= limit; k++) {
            long kk = k;
            rem[k] = kk * kk - kk + 1;
        }

        for (int p : primes) {
            if (p == 2)
                continue;
            List<Integer> roots = new ArrayList<>();
            if (p == 3) {
                roots.add(2);
            } else {
                long mod = p;
                long negThree = (mod + mod - 3) % mod;
                if (powMod(negThree, (mod - 1) / 2, mod) != 1)
                    continue;

                long sq = tonelliShanks(negThree, mod);
                long inv2 = (mod + 1) / 2;
                int r1 = (int) mulMod((1 + sq) % mod, inv2, mod);
                int r2 = (int) mulMod((1 + mod - sq) % mod, inv2, mod);
                roots.add(r1);
                if (r2 != r1) {
                    roots.add(r2);
                }
            }

            for (int root : roots) {
                int k = root;
                if (k == 0)
                    k += p;
                for (; k <= limit; k += p) {
                    long value = rem[k];
                    if (value % p != 0)
                        continue;
                    while (value % p == 0) {
                        value /= p;
                    }
                    rem[k] = value;
                    lpfQ[k] = p;
                }
            }
        }

        long total = 0;
        for (int k = 1; k <= limit; k++) {
            long remVal = rem[k];
            if (remVal > 1) {
                lpfQ[k] = Math.max(lpfQ[k], remVal);
            }
            long lpfKp1 = lpfSmall[k + 1];
            long largest = Math.max(lpfKp1, lpfQ[k]);
            total += (largest - 1);
        }
        return total;
    }

    public static String solve() {
        long ans = solve(2000000);
        return String.valueOf(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}