Problem 245: Coresilience

Project Euler Problem 245 Solution

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

Problem Summary Define the coresilience of a positive integer by $$C(n)=\frac{n-\varphi(n)}{n-1}.$$ Problem 245 asks for the sum of all composite numbers $n\le 2\cdot 10^{11}$ whose coresilience is a unit fraction. Equivalently, we must find every composite $n$ for which $$\frac{n-1}{n-\varphi(n)}\in\mathbb Z.$$ If we write $$n-1=k(n-\varphi(n)),\qquad k\in\mathbb Z_{\ge 2},$$ then $C(n)=1/k$. A direct sweep up to $2\cdot 10^{11}$ with explicit totients would be far too slow, so the solution turns this divisibility condition into a structural classification of the admissible composite numbers. Mathematical Approach The integer $k$ is the natural organizing parameter of the search. Rearranging the defining equation gives the invariant $$k\varphi(n)=(k-1)n+1.$$ Everything in the final algorithm comes from forcing this identity to coexist with the prime factorization of $n$. Squarefree structure is forced Suppose $p^2\mid n$ for some prime $p$. Then $p\mid n$, and also $p\mid \varphi(n)$, because the totient of any integer divisible by $p^2$ still carries a factor of $p$. Hence $p\mid n-\varphi(n)$. But $n\equiv 0\pmod p$ implies $n-1\equiv -1\pmod p$, so $p\nmid n-1$. Therefore $n-\varphi(n)$ cannot divide $n-1$....

Detailed mathematical approach

Problem Summary

Euler's totient $\varphi(n)$ counts the integers $1\le j\le n$ that have no prime factor in common with $n$. A fraction that cannot be cancelled down is called resilient, and the resilience $R(d)$ of a denominator $d\gt1$ is the share of its proper fractions $1/d,\dots,(d-1)/d$ that are resilient, $R(d)=\varphi(d)/(d-1)$; for example $R(12)=4/11$. The coresilience of a number $n\gt1$ is

$$C(n)=\frac{n-\varphi(n)}{n-1}.$$

Its numerator $n-\varphi(n)$ counts the integers $1\le j\le n$ that share a prime factor with $n$. For a prime $p$ only $j=p$ does, so $C(p)=1/(p-1)$ is always a unit fraction, a fraction with numerator $1$. The task is to find the sum of all composite integers $n$ with $1\lt n\le2\times10^{11}$ for which $C(n)$ is a unit fraction. ^[1] We write $L=2\times10^{11}$ for the limit.

The whole solution in one sentence: writing $C(n)=1/k$, every solution is a product of distinct primes larger than $k$; with two primes $p,q$ the condition becomes $(p-k)(q-k)=k^2-k+1$, so these solutions come from the divisors of $k^2-k+1$, and with three or more primes the largest prime is fixed by the others, so a search over the smaller primes finds the rest.

Mathematical Approach

Computing $\varphi(n)$ for all $2\times10^{11}$ numbers is out of reach, so we first find out what the prime factorisation of a solution must look like; then the solutions can be enumerated directly.

1. Small cases

The first composite numbers give $C(4)=2/3$, $C(6)=4/5$, $C(8)=4/7$, $C(9)=3/8$, $C(10)=2/3$, $C(12)=8/11$ and $C(14)=8/13$, none of them a unit fraction. The first unit fraction is $C(15)=(15-8)/14=1/2$. A direct search up to $1000$ finds exactly six composite numbers with this property; they begin the sequence A160599 of the On-Line Encyclopedia of Integer Sequences: ^[2]

$15=3\cdot5$ with $C(15)=1/2$;
$85=5\cdot17$ with $C(85)=1/4$;
$255=3\cdot5\cdot17$ with $C(255)=1/2$;
$259=7\cdot37$ with $C(259)=1/6$;
$391=17\cdot23$ with $C(391)=1/10$;
$589=19\cdot31$ with $C(589)=1/12$.

Their sum is $1594$. Four patterns stand out: every number is odd, no prime factor appears twice, every prime factor is larger than the denominator $k$ of $C(n)=1/k$, and that denominator is even. The next sections prove all four.

2. The integer k

Write $C(n)=1/k$. Then $n-1=k\,(n-\varphi(n))$, and rearranging gives the equation behind everything below:

$$k\,\varphi(n)=(k-1)\,n+1.\qquad(1)$$

Here $k\ge2$, because $k=1$ would mean $\varphi(n)=1$, which happens only for $n\le2$. For example, $\varphi(85)=64$ and $4\cdot64=256=3\cdot85+1$, so $k=4$ for $n=85$. We use two facts about the totient: ^[3] $\varphi(p^e)=p^{e-1}(p-1)$ for a prime power, and $\varphi(ab)=\varphi(a)\varphi(b)$ for coprime $a,b$; so $\varphi(n)=\prod(p_i-1)$ for a product $n=p_1\cdots p_r$ of distinct primes.

3. No prime factor appears twice

Suppose that $p^2$ divides $n$. Then $p$ divides $\varphi(n)$, which contains the factor $p^{e-1}$ with $e\ge2$, and $p$ divides $n$, so $p$ divides $n-\varphi(n)$. That number divides $n-1$, so $p$ would divide both $n$ and $n-1$, which is impossible. Hence a solution is squarefree: $n=p_1p_2\cdots p_r$ with $r\ge2$ distinct primes. For example $n=9$ has $n-\varphi(n)=3$, which does not divide $n-1=8$. In Lehmer's totient problem, a close relative, no composite $n$ with $\varphi(n)\mid n-1$ is known, and Lehmer showed that one would have to be odd and squarefree. ^[4]

4. Every prime factor exceeds k

Dividing (1) by $kn$ gives

$$\frac{\varphi(n)}{n}=1-\frac1k+\frac{1}{kn}\gt1-\frac1k.$$

For squarefree $n$ the left side is the product of the numbers $1-1/p$ over the primes $p\mid n$. Each of them is below $1$, so each single factor is at least as large as the whole product: $1-1/p\gt1-1/k$, that is, $p\gt k$ for every prime factor $p$. Three consequences follow. Since $k\ge2$, every prime factor is odd, so $n$ is odd. Then $\varphi(n)=\prod(p_i-1)$ is even, so $n-\varphi(n)$ is odd while $n-1$ is even, and $k=(n-1)/(n-\varphi(n))$ must be even. Finally, the size of $n$ bounds $k$: two prime factors above $k$ give $n\gt k^2$, and three or more give $n\ge(k+1)(k+2)(k+3)$, which exceeds $(k+1)^3$.

5. Two primes: a divisor problem

Let $n=pq$ with primes $p\lt q$. Then $\varphi(n)=(p-1)(q-1)$ and $n-\varphi(n)=p+q-1$, so the condition reads $pq-1=k(p+q-1)$. Rearranging and adding $k^2$ to both sides gives

$$(p-k)(q-k)=k^2-k+1.$$

We write $K(k)=k^2-k+1$. So $a=p-k$ and $b=q-k$ are complementary divisors, $ab=K(k)$, and conversely every divisor pair for which $k+a$ and $k+b$ are both prime gives a solution, because the algebra can be read backwards. Since $(k-1)^2\lt K(k)\lt k^2$, $K(k)$ is never a square, so $a\ne b$, and we take $a\lt b$. The even values of $k$ up to $12$:

$k=2$: $K=3$, the pair $(1,3)$ gives $3\cdot5=15$;
$k=4$: $K=13$, the pair $(1,13)$ gives $5\cdot17=85$;
$k=6$: $K=31$, the pair $(1,31)$ gives $7\cdot37=259$;
$k=8$: $K=57=3\cdot19$, the pairs $(1,57)$ and $(3,19)$ give $9\cdot65$ and $11\cdot27$, and neither is a pair of primes;
$k=10$: $K=91=7\cdot13$, the pairs $(1,91)$ and $(7,13)$ give $11\cdot101=1111$ and $17\cdot23=391$;
$k=12$: $K=133=7\cdot19$, the pair $(1,133)$ gives $13\cdot145$ with $145=5\cdot29$, and $(7,19)$ gives $19\cdot31=589$.

For odd $k$ every divisor of the odd number $K(k)$ is odd, so $k+a$ is even and larger than $2$: no solution, as Section 4 predicts. The condition also limits $k$: it gives $n=k(p+q)-k+1$, and $p+q=2k+a+b$ with $a+b\ge2\sqrt{K(k)}\gt2k-2$, so

$$\begin{aligned}n&\gt k(4k-2)-k+1\\&=4k^2-3k+1\gt(2k-1)^2.\end{aligned}$$

For $L=2\times10^{11}$ only $k\le223\,607$ can contribute. The program uses the simpler bound $n\gt k^2$ and tries every $k\le\sqrt L$.

6. Which primes divide $k^2-k+1$

To factorise all $K(k)$ at once, we need to know which $k$ make a given prime $p$ a divisor. $K(k)=k(k-1)+1$ is odd, so $p=2$ never divides it, and multiplying by $4$ completes the square:

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

For $p=3$ this shows that $3\mid K(k)$ exactly when $3\mid2k-1$, that is, $k\equiv2\pmod3$. For a prime $p\gt3$, $p\mid K(k)$ holds exactly when $(2k-1)^2\equiv-3\pmod p$. By quadratic reciprocity, $-3$ is a square modulo $p$ exactly when $p\equiv1\pmod3$, ^[5] so the primes $p\equiv2\pmod3$ never divide any $K(k)$. For $p\equiv1\pmod3$, let $s$ be a square root of $-3$ modulo $p$. Then $2k-1\equiv\pm s$, and $p\mid K(k)$ holds for exactly the two residue classes

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

which differ because $s\not\equiv0$; dividing by $2$ means multiplying by $(p+1)/2$, the inverse of $2$ modulo $p$. For $p=7$: $-3\equiv4=2^2$, so $s=2$, and $1/2\equiv4$ gives $k\equiv3$ and $k\equiv5\pmod7$; indeed $K(3)=7$, $K(5)=21$, $K(10)=91$ and $K(12)=133$ are multiples of $7$. For $p=13$: $6^2\equiv-3$ and $1/2\equiv7$ give $k\equiv4$ and $k\equiv10\pmod{13}$, matching $K(4)=13$ and $K(10)=91$. For large primes the square root comes from the Tonelli–Shanks algorithm, which builds it from a few modular powers; for $p\equiv3\pmod4$ it is simply $(-3)^{(p+1)/4}\bmod p$. ^[6]

7. Factorising all values with a sieve

The sieve of Eratosthenes crosses out the multiples of one prime after another. ^[7] In the same spirit, the program keeps for every $k\le\sqrt L$ a remainder that starts at $K(k)$, and for every prime $p\le\sqrt L$ that can divide some $K(k)$ it walks only through the $k$ of its residue classes, $k=r,r+p,r+2p,\dots$, divides every factor $p$ out of the remainder and records the exponent, instead of trying every prime on every value. When all primes up to $\sqrt L$ have been used, each remainder is $1$ or a single prime: two leftover prime factors would both exceed $\sqrt L\ge k$, and their product would exceed $k^2\gt K(k)$. For $k=10$ the walks of $7$ and $13$ reduce $91$ to $1$; for the real limit, $306\,847$ of the $447\,212$ values keep such a large prime. From the factorisation $K=\prod p_i^{e_i}$ the divisors follow one prime at a time, by multiplying the divisors found so far by $p_i,p_i^2,\dots,p_i^{e_i}$.

8. Three or more primes: the last prime is forced

Now let $n$ have $r\ge3$ prime factors. Write $n=mp$, where $p$ is the largest prime factor and $m=q_1q_2\cdots q_t$ the product of the others, $t=r-1\ge2$, and let $A=\varphi(m)=(q_1-1)\cdots(q_t-1)$. Then $\varphi(n)=A(p-1)$, and (1) becomes $kA(p-1)=(k-1)mp+1$. Collecting the terms with $p$ gives

$$\begin{aligned}p\,D&=kA+1,\\D&=kA-(k-1)\,m.\end{aligned}$$

Once $k$ and the smaller primes are chosen, the largest prime is therefore determined: $D$ must be positive and divide $kA+1$, and $p=(kA+1)/D$ must be a prime larger than $q_t$. Two examples:

$k=18$, $m=43\cdot47=2021$, $A=42\cdot46=1932$: here $kA+1=34777$ and $D=34776-34357=419$, so $p=34777/419=83$, a prime. Indeed $n=43\cdot47\cdot83=167743$ has $\varphi(n)=158424$, $n-\varphi(n)=9319$ and $18\cdot9319=n-1$.
$k=2$, $m=15$, $A=8$: here $D=16-15=1$ and $p=17$, which gives $255$.

The second example is a special case. If $m$ itself satisfies (1) with the same $k$, that is, $kA=(k-1)m+1$, then $D=1$ and $p=kA+1$, and $mp$ satisfies (1) again whenever $p$ is prime, because $\varphi(mp)=A\cdot kA$ and

$$\begin{aligned}&(k-1)\,m\,(kA+1)+1\\&\quad=kA\bigl((k-1)m+1\bigr)=k^2A^2.\end{aligned}$$

For $k=2$ this gives the chain $15\to255\to65535\to4294967295$: the new factors $17$, $257$ and $65537$ are the Fermat primes $F_j=2^{2^j}+1$, and the members are $F_0F_1\cdots F_j=F_{j+1}-2$. The chain stops because the next factor $2^{32}+1=F_5=641\cdot6700417$ is not prime. ^[8] The OEIS entry notes that every product of two or more consecutive Fermat primes belongs to the sequence; ^[2] all ten such products of the five known Fermat primes lie below $L$.

9. Searching for the smaller primes

For each $k$ we build $m$ one prime at a time in a depth-first search: a path $q_1\lt q_2\lt\cdots$ of primes above $k$, with $m$ and $A$ carried along, and at every node with $t\ge2$ we test the forced prime of Section 8. Two size bounds limit the search, and a sign rule can prune it further.

Only $k$ with $(k+1)^3\lt L$ can occur, by Section 4; for $L=2\times10^{11}$ these are $k\le5847$.
A prime $q$ is appended to $m$ only if $mq^2\le L$, because the final number contains $q$ and at least one larger prime. In particular all $q_i$ are at most $\sqrt L$.
A path whose $D$ is not positive can be abandoned, for the reason given next; the Java program does so, the C++ program does not (see How the Code Works).

Appending a prime $q$ to $m$ changes the denominator to

$$\begin{aligned}D(mq)&=kA(q-1)-(k-1)\,mq\\&=q\,D(m)-kA,\end{aligned}$$

so once $D\le0$, every extension of the path has $D\lt0$ and can never produce a prime. The condition $D\gt0$ says $\varphi(m)/m\gt(k-1)/k$, the inequality of Section 4 for the part $m$. With two primes, $D=(q_1-k)(q_2-k)-k(k-1)$. For $k=18$ the start $19,23$ gives $D=1\cdot5-306=-301$, so no number that begins with $19\cdot23$ can work, while the start $43,47$ gives the $D=419$ of Section 8.

10. Every solution is found exactly once

Take any solution $n$. Its $k$ is determined by $n$, and its prime factors are distinct by Section 3. With two primes, $n\gt k^2$ gives $k\le\sqrt L$, and $(p-k,q-k)$ is a divisor pair of $K(k)$ from Section 7. With $r\ge3$ primes, $k\le5847$, the path of the $r-1$ smallest primes obeys $mq^2\le L$ at every step, and every part of it has $D\gt0$, because the complete path has $D=(kA+1)/p\gt0$ and a non-positive $D$ never becomes positive again; so the search reaches the path, and Section 8 returns the largest prime. Conversely, every recorded number is a product of distinct primes that satisfies (1). Since $k$, the divisor pair and the path are all determined by $n$, no number is recorded twice.

How the Code Works

The C++ program ^[9] accepts an optional limit after -n, with $2\times10^{11}$ as the default, and checks itself first unless --no-validate is given. solve_fast computes $\lfloor\sqrt L\rfloor=447\,213$ with isqrt_u64, and sieve_primes lists the primes up to that bound with the sieve of Eratosthenes, keeping the sieve as a lookup table. It then calls collect_semiprimes and collect_multi_primes, sorts the collected numbers, removes duplicates as a safeguard and adds them up in an unsigned 128-bit integer.

collect_semiprimes implements Sections 5 to 7. The vector rem starts with $k^2-k+1$ for $2\le k\le447\,213$, and factors[k] collects pairs of a prime and its exponent. The prime $2$ and the primes $p\equiv2\pmod3$ are skipped; for $p=3$ the loop visits $k=2,5,8,\dots$; for $p\equiv1\pmod3$ it calls tonelli_shanks(p - 3, p) for a square root of $-3$, forms the two classes with the inverse $(p+1)/2$ of $2$, and sieve_root walks each class and divides the factors $p$ out. A remainder above $1$ is then recorded as a prime with exponent $1$. For each $k$ the divisors are built from the factor list, and every divisor $a\le b=K/a$ gives $p=k+a$ and $q=k+b$; the pair is skipped if $p=q$ or if the 128-bit product $pq$ exceeds $L$, and $pq$ is recorded when both numbers are prime.

collect_multi_primes runs through $k=2,\dots,\lfloor L^{1/3}\rfloor-1$, with $\lfloor L^{1/3}\rfloor=5848$ from icbrt_u64, and starts dfs_multi with $m=A=1$ at the first prime above $k$. Each call computes $D$ and $kA+1$; if $D$ divides $kA+1$, the quotient $p$ is recorded as $mp$ when it exceeds the last prime, $mp\le L$, at least two primes have been chosen and $p$ is prime. The call then tries the following primes $q$ in increasing order, stops at the first with $mq^2\gt L$, and recurses with $mq$, $A(q-1)$ and $q$. The program computes $D$ in unsigned 128-bit arithmetic, so its test denom <= 0 fires only for $D=0$, which is impossible: $kA=(k-1)m$ would make $k$ divide $m$, since $k$ and $k-1$ are coprime, but every prime factor of $m$ exceeds $k$. A negative $D$ wraps around to a number close to $2^{128}$, and the remainder of $kA+1\lt10^{16}$ is then $kA+1$ itself, so no wrong number is recorded; but the dead branches of Section 9 are still searched, with $1\,481\,541\,274$ calls instead of the $1\,393\,011\,937$ of a signed test, about $6\%$ more.

is_prime answers from the sieve table below $447\,214$ and otherwise calls is_prime_mr, which divides by the primes up to $37$ and then runs the Miller–Rabin test: with $n-1=2^sd$ and $d$ odd, a base $a$ passes if $a^d\equiv1$ or $a^{2^jd}\equiv-1\pmod n$ for some $j\lt s$, and $n$ counts as prime when the bases $2,3,5,7,11,13$ all pass. mod_mul forms products modulo $n$ in 128 bits and mod_pow computes powers by repeated squaring. These six bases are known to give the correct answer for every number below $3\,474\,749\,660\,383$, ^[10] and the program only tests numbers up to $L$. run_validation compares solve_fast with brute_sum, which builds a table of all totients with a sieve and adds every non-prime $n\ge4$ for which $n-\varphi(n)$ divides $n-1$.

The Java program follows the same plan in one method solve, without checkpoints; it tests products by division, as in pp > LIMIT / qq, collects the results in a HashSet and uses BigInteger for the modular arithmetic of Tonelli–Shanks and Miller–Rabin, with the same six bases. Its dfsMulti computes $D$ as a signed long, so there the test $D\le0$ does cut off dead branches. Its factor table stores primes as int, so the cast (int) rem[k] garbles a leftover prime above $2^{31}-1$, as happens for $82\,586$ values of $k$. No solution is lost, since a divisor $a\lt\sqrt K\lt k$ never contains the leftover prime and $b=K/a$ comes from the exact $K$; the garbled factor only adds candidates $a$ that do not divide $K$, and a simulation of the Java arithmetic shows that none of them passes the tests at this limit. The Python file is a thin bridge: it compiles solutionsCpp/Euler245.cpp with clang++ or g++ when its hidden binary is missing or older than the source, runs it with the limit $2\times10^{11}$ and --no-validate, and returns the last line of the output.

Complexity and Verification

The prime sieve up to $\sqrt L$ takes $O(\sqrt L\log\log L)$ time and finds $37\,499$ primes. The residue-class sieve of Section 7 visits $977\,402$ values of $k$ in all; this is $O(\sqrt L\log\log L)$, because a prime $p$ has at most $2\sqrt L/p+2$ visits and the sum of $1/p$ over the primes up to $x$ grows like $\log\log x$. On top of that come $18\,718$ calls of Tonelli–Shanks, one for each prime $p\equiv1\pmod3$ below $\sqrt L$, each costing a few modular powers. The $447\,212$ values $K(k)$ have $4\,634\,348$ divisors in total, at most $288$ for one value, that is $2\,317\,174$ pairs, of which $256\,044$ satisfy $pq\le L$ and reach the primality tests. This whole stage takes about a tenth of a second.

The search of Section 9 is the expensive part. Besides the $5846$ initial calls, the C++ program makes $216\,814\,048$ calls of dfs_multi at depth one, $988\,235\,862$ at depth two and $276\,485\,518$ deeper. Of the depth-one calls, $214\,735\,272$ end at once: the outer loop tries every prime up to $\sqrt L$ as $q_1$, although a smallest prime above $L^{1/3}$ leaves no room for two larger ones. The depth-two calls, one for every $k$ and every pair $q_1\lt q_2$ of primes above $k$ with $q_1q_2^2\le L$, are the real work; their number grows roughly like $L/(\log L)^2$, from $7\,756\,219$ for $L=10^9$ to $988\,235\,862$ for $L=2\times10^{11}$, $127$ times as many for a $200$ times larger limit. About half of all calls belong to odd $k$, which by Section 4 can never succeed; the program does not use this. Each call performs a 128-bit division; the search takes about eight seconds on a desktop computer, nearly all of the running time. The memory is $O(\sqrt L\log\log L)$ words for the remainders, the sieve and the factor lists, whose total length is the number of sieve visits plus at most one leftover prime per value; in practice the C++ program needs about 50 MB.

Before the main computation, the C++ program compares brute_sum and solve_fast with $1594$, $595394$ and $11149065$ for $L=10^3$, $10^5$ and $10^6$; the first value is the sum of the six numbers of Section 1. For this article, a separate brute force with a totient sieve found $547$ qualifying numbers up to $10^9$ with the sum $145645009845$, which is also what the C++ program returns for that limit. An independent program that factorises every $K(k)$ by trial division, without the residue classes of Section 6, and runs its own search with signed arithmetic finds the same $4877$ numbers up to $2\times10^{11}$, namely $4830$ with two prime factors, $38$ with three, $7$ with four and $2$ with five, and their sum agrees with the value printed by all three published programs. Instrumented copies of the C++ program confirm the claims above: no number is found twice, every leftover factor is prime, and the largest $k$ that occurs is $221\,940$ for two prime factors and $1550$ for more. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The reduction to divisor pairs of $k^2-k+1$, the forced last prime, the pruning rules of the search and the operation counts are derived explicitly in this article.

Project Euler 245 — Coresilience. The official statement defines the resilience R(d) = φ(d)/(d − 1), with R(12) = 4/11, and the coresilience C(n) = (n − φ(n))/(n − 1), notes that C(p) = 1/(p − 1) for a prime p, and asks for the sum of all composite n with 1 < n ≤ 2 × 10¹¹ for which C(n) is a unit fraction.
OEIS A160599 — Composite numbers k for which k − phi(k) divides k − 1. The entry lists 15, 85, 255, 259, 391, 589, 1111, … and remarks that the sequence contains every product of successive Fermat primes.
Wikipedia — Euler's totient function. The article gives Euler's product φ(n) = n ∏(1 − 1/p) over the distinct primes p dividing n, the value φ(pᵏ) = pᵏ⁻¹(p − 1) for prime powers, and φ(mn) = φ(m)φ(n) for coprime m and n.
Wikipedia — Lehmer's totient problem. The unsolved question whether some composite n has φ(n) dividing n − 1; Lehmer showed that such an n must be odd and squarefree with at least seven distinct prime factors.
Wikipedia — Quadratic reciprocity. For primes p > 3, the section on q = ±3 shows that −3 is a quadratic residue modulo p exactly when p ≡ 1 (mod 3); Section 6 applies this to 4K(k) = (2k − 1)² + 3 and treats 2 and 3 separately.
Wikipedia — Tonelli–Shanks algorithm. The algorithm solves r² ≡ n (mod p) for an odd prime p and a quadratic residue n with the help of a quadratic non-residue; for p ≡ 3 (mod 4) the solutions are ±n^((p+1)/4).
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime as composite, starting from the square of the prime, and runs in O(n log log n) time; Section 7 applies the same idea to the values k² − k + 1.
Wikipedia — Fermat number. Fₙ = 2^(2ⁿ) + 1; F₀, …, F₄ = 3, 5, 17, 257, 65537 are the only known Fermat primes, Fₙ = F₀F₁⋯Fₙ₋₁ + 2 for n ≥ 1, and Euler found F₅ = 641 × 6700417.
C++ — Euler245.cpp. The linked immutable C++ revision contains collect_semiprimes, dfs_multi, collect_multi_primes, solve_fast and the checks in run_validation. It is the implementation record for this article.
Wikipedia — Miller–Rabin primality test. With n − 1 = 2ˢ·d and d odd, the test checks whether n is a strong probable prime to a base a; the article lists that testing the bases 2, 3, 5, 7, 11 and 13 is enough for every n < 3,474,749,660,383.

Problem 245 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = uint64_t;
using u128 = __uint128_t;

constexpr u64 kDefaultLimit = 200000000000ULL;

u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1) * (r + 1) <= x) ++r;
    while (r * r > x) --r;
    return r;
}

u64 icbrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
    while ((r + 1) * (r + 1) * (r + 1) <= x) ++r;
    while (r * r * r > x) --r;
    return r;
}

std::vector<int> sieve_primes(int max_n, std::vector<bool>* is_prime_out) {
    std::vector<bool> is_prime(max_n + 1, true);
    is_prime[0] = false;
    if (max_n >= 1) is_prime[1] = false;
    for (int i = 2; (int64_t)i * i <= max_n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= max_n; j += i) is_prime[j] = false;
    }
    std::vector<int> primes;
    for (int i = 2; i <= max_n; ++i) {
        if (is_prime[i]) primes.push_back(i);
    }
    if (is_prime_out) *is_prime_out = std::move(is_prime);
    return primes;
}

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

u64 mod_pow(u64 a, u64 e, u64 mod) {
    u64 result = 1 % mod;
    u64 base = a % mod;
    while (e > 0) {
        if (e & 1) result = mod_mul(result, base, mod);
        base = mod_mul(base, base, mod);
        e >>= 1;
    }
    return result;
}

bool is_prime_mr(u64 n) {
    if (n < 2) return false;
    static const u64 small_primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    for (u64 p : small_primes) {
        if (n % p == 0) return n == p;
    }
    u64 d = n - 1;
    int s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        ++s;
    }
    auto witness = [&](u64 a) {
        if (a % n == 0) return false;
        u64 x = mod_pow(a, d, n);
        if (x == 1 || x == n - 1) return false;
        for (int i = 1; i < s; ++i) {
            x = mod_mul(x, x, n);
            if (x == n - 1) return false;
        }
        return true;
    };
    static const u64 bases[] = {2, 3, 5, 7, 11, 13};
    for (u64 a : bases) {
        if (witness(a)) return false;
    }
    return true;
}

bool is_prime(u64 n, const std::vector<bool>& is_prime_small) {
    if (n < is_prime_small.size()) return is_prime_small[n];
    return is_prime_mr(n);
}

u64 tonelli_shanks(u64 n, u64 p) {
    if (n == 0) return 0;
    if (p == 2) return 1;
    if (p % 4 == 3) return mod_pow(n, (p + 1) / 4, p);
    u64 q = p - 1;
    int s = 0;
    while ((q & 1) == 0) {
        q >>= 1;
        ++s;
    }
    u64 z = 2;
    while (mod_pow(z, (p - 1) / 2, p) != p - 1) ++z;
    u64 c = mod_pow(z, q, p);
    u64 x = mod_pow(n, (q + 1) / 2, p);
    u64 t = mod_pow(n, q, p);
    int m = s;
    while (t != 1) {
        int i = 1;
        u64 t2i = mod_mul(t, t, p);
        while (t2i != 1) {
            t2i = mod_mul(t2i, t2i, p);
            ++i;
        }
        u64 b = mod_pow(c, 1ULL << (m - i - 1), p);
        x = mod_mul(x, b, p);
        u64 b2 = mod_mul(b, b, p);
        t = mod_mul(t, b2, p);
        c = b2;
        m = i;
    }
    return x;
}

struct Factor {
    u64 prime = 0;
    int exp = 0;
};

void collect_semiprimes(u64 limit, int threads, const std::vector<int>& primes,
                        const std::vector<bool>& is_prime_small, std::vector<u64>& out) {
    u64 k_max = isqrt_u64(limit);
    std::vector<u64> rem(k_max + 1);
    std::vector<std::vector<Factor>> factors(k_max + 1);
    for (u64 k = 2; k <= k_max; ++k) {
        rem[k] = k * k - k + 1;
    }

    for (int p : primes) {
        if (p > static_cast<int>(k_max)) break;
        if (p == 2) continue;
        if (p == 3) {
            for (u64 k = 2; k <= k_max; k += 3) {
                if (rem[k] % 3 != 0) continue;
                int exp = 0;
                while (rem[k] % 3 == 0) {
                    rem[k] /= 3;
                    ++exp;
                }
                factors[k].push_back({3, exp});
            }
            continue;
        }
        if (p % 3 != 1) continue;

        u64 n = static_cast<u64>(p - 3);
        u64 root = tonelli_shanks(n, static_cast<u64>(p));
        u64 inv2 = (static_cast<u64>(p) + 1) / 2;
        u64 r1 = (1 + root) % p;
        r1 = (r1 * inv2) % p;
        u64 r2 = (1 + p - root) % p;
        r2 = (r2 * inv2) % p;

        auto sieve_root = [&](u64 r) {
            if (r < 2) {
                u64 add = (2 - r + p - 1) / p;
                r += add * static_cast<u64>(p);
            }
            for (u64 k = r; k <= k_max; k += p) {
                if (rem[k] % p != 0) continue;
                int exp = 0;
                while (rem[k] % p == 0) {
                    rem[k] /= p;
                    ++exp;
                }
                factors[k].push_back({static_cast<u64>(p), exp});
            }
        };

        sieve_root(r1);
        if (r2 != r1) sieve_root(r2);
    }

    for (u64 k = 2; k <= k_max; ++k) {
        if (rem[k] > 1) factors[k].push_back({rem[k], 1});
    }

    int use_threads = std::max(1, threads);
    std::vector<std::vector<u64>> locals(use_threads);
    auto worker = [&](int tid, u64 start, u64 end) {
        std::vector<u64> divisors;
        divisors.reserve(64);
        for (u64 k = start; k < end; ++k) {
            u64 K = k * k - k + 1;
            divisors.clear();
            divisors.push_back(1);
            for (const auto& f : factors[k]) {
                size_t size = divisors.size();
                u64 mult = 1;
                for (int e = 1; e <= f.exp; ++e) {
                    mult *= f.prime;
                    for (size_t i = 0; i < size; ++i) {
                        divisors.push_back(divisors[i] * mult);
                    }
                }
            }
            for (u64 a : divisors) {
                u64 b = K / a;
                if (a > b) continue;
                u64 p = k + a;
                u64 q = k + b;
                if (p == q) continue;
                if (static_cast<u128>(p) * q > limit) continue;
                if (is_prime(p, is_prime_small) && is_prime(q, is_prime_small)) {
                    locals[tid].push_back(p * q);
                }
            }
        }
    };

    std::vector<std::thread> pool;
    pool.reserve(use_threads);
    u64 total = k_max - 1;
    u64 chunk = (total + use_threads - 1) / use_threads;
    for (int t = 0; t < use_threads; ++t) {
        u64 start = 2 + t * chunk;
        u64 end = std::min<u64>(2 + (t + 1) * chunk, k_max + 1);
        if (start >= end) break;
        pool.emplace_back(worker, t, start, end);
    }
    for (auto& th : pool) th.join();

    for (auto& vec : locals) {
        out.insert(out.end(), vec.begin(), vec.end());
    }
}

void dfs_multi(int k, size_t idx, u64 m, u64 A, int last_prime, int depth, u64 limit,
               const std::vector<int>& primes, const std::vector<bool>& is_prime_small,
               std::vector<u64>& out) {
    u128 denom = static_cast<u128>(k) * A - static_cast<u128>(k - 1) * m;
    if (denom <= 0) return;
    u128 numer = static_cast<u128>(k) * A + 1;
    if (numer % denom == 0) {
        u64 p = static_cast<u64>(numer / denom);
        if (p > static_cast<u64>(last_prime) && static_cast<u128>(m) * p <= limit) {
            if (depth >= 2 && is_prime(p, is_prime_small)) {
                out.push_back(m * p);
            }
        }
    }

    for (size_t i = idx; i < primes.size(); ++i) {
        int q = primes[i];
        if (q <= k) continue;
        if (static_cast<u128>(m) * q * q > limit) break;
        dfs_multi(k, i + 1, m * static_cast<u64>(q), A * static_cast<u64>(q - 1), q,
                  depth + 1, limit, primes, is_prime_small, out);
    }
}

void collect_multi_primes(u64 limit, int threads, const std::vector<int>& primes,
                          const std::vector<bool>& is_prime_small, std::vector<u64>& out) {
    u64 k_root = icbrt_u64(limit);
    if (k_root == 0) return;
    int k_max = static_cast<int>(k_root) - 1;
    if (k_max < 2) return;

    int use_threads = std::max(1, threads);
    std::vector<std::vector<u64>> locals(use_threads);
    auto worker = [&](int tid, int start_k, int end_k) {
        for (int k = start_k; k < end_k; ++k) {
            auto it = std::upper_bound(primes.begin(), primes.end(), k);
            size_t idx = static_cast<size_t>(it - primes.begin());
            dfs_multi(k, idx, 1, 1, 1, 0, limit, primes, is_prime_small, locals[tid]);
        }
    };

    int total = k_max - 1;
    int chunk = (total + use_threads - 1) / use_threads;
    std::vector<std::thread> pool;
    for (int t = 0; t < use_threads; ++t) {
        int start_k = 2 + t * chunk;
        int end_k = std::min(2 + (t + 1) * chunk, k_max + 1);
        if (start_k >= end_k) break;
        pool.emplace_back(worker, t, start_k, end_k);
    }
    for (auto& th : pool) th.join();

    for (auto& vec : locals) {
        out.insert(out.end(), vec.begin(), vec.end());
    }
}

u64 solve_fast(u64 limit, int threads) {
    u64 k_max = isqrt_u64(limit);
    std::vector<bool> is_prime_small;
    std::vector<int> primes = sieve_primes(static_cast<int>(k_max), &is_prime_small);
    std::vector<u64> results;
    results.reserve(1024);
    collect_semiprimes(limit, threads, primes, is_prime_small, results);
    collect_multi_primes(limit, threads, primes, is_prime_small, results);
    std::sort(results.begin(), results.end());
    results.erase(std::unique(results.begin(), results.end()), results.end());
    u128 sum = 0;
    for (u64 n : results) sum += n;
    return static_cast<u64>(sum);
}

u64 brute_sum(u64 limit) {
    std::vector<u64> phi(limit + 1);
    for (u64 i = 0; i <= limit; ++i) phi[i] = i;
    for (u64 i = 2; i <= limit; ++i) {
        if (phi[i] != i) continue;
        for (u64 j = i; j <= limit; j += i) {
            phi[j] -= phi[j] / i;
        }
    }
    u64 sum = 0;
    for (u64 n = 4; n <= limit; ++n) {
        u64 c = n - phi[n];
        if ((n - 1) % c != 0) continue;
        if (phi[n] == n - 1) continue;
        sum += n;
    }
    return sum;
}

void run_validation() {
    struct Check {
        u64 limit;
        u64 expected;
    };
    const Check checks[] = {
        {1000, 1594},
        {100000, 595394},
        {1000000, 11149065},
    };
    for (const auto& chk : checks) {
        u64 brute = brute_sum(chk.limit);
        if (brute != chk.expected) {
            std::cerr << "Validation failed (brute) for N=" << chk.limit
                      << ": got " << brute << " expected " << chk.expected << "\n";
            std::exit(1);
        }
        u64 fast = solve_fast(chk.limit, 1);
        if (fast != chk.expected) {
            std::cerr << "Validation failed (fast) for N=" << chk.limit
                      << ": got " << fast << " expected " << chk.expected << "\n";
            std::exit(1);
        }
    }
    std::cerr << "Validation checkpoints passed.\n";
}

void usage(const char* argv0) {
    std::cerr << "Usage: " << argv0 << " [-n LIMIT] [-t THREADS] [--no-validate]\n";
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    u64 limit = kDefaultLimit;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;
    bool validate = true;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "-n" && i + 1 < argc) {
            limit = std::stoull(argv[++i]);
        } else if (arg == "-t" && i + 1 < argc) {
            threads = std::max(1, std::stoi(argv[++i]));
        } else if (arg == "--no-validate") {
            validate = false;
        } else if (arg == "--validate") {
            validate = true;
        } else {
            usage(argv[0]);
            return 1;
        }
    }

    if (validate) run_validation();

    u64 answer = solve_fast(limit, threads);
    std::cout << answer << "\n";
    return 0;
}

Python

from __future__ import annotations

import os
import shutil
import subprocess
from pathlib import Path


def pick_compiler() -> str:
    for compiler in ("clang++", "g++"):
        if shutil.which(compiler):
            return compiler
    raise RuntimeError("No C++ compiler found (clang++/g++).")


def ensure_bridge_binary(root: Path) -> Path:
    src = root / "solutionsCpp" / "Euler245.cpp"
    bin_path = root / "solutionsCpp" / ".euler245_py_bridge"
    if (not bin_path.exists()) or (src.stat().st_mtime > bin_path.stat().st_mtime):
        compiler = pick_compiler()
        cmd = [compiler, "-std=c++17", "-O2", str(src), "-o", str(bin_path)]
        subprocess.check_call(cmd)
    return bin_path


def solve() -> str:
    root = Path(__file__).resolve().parent.parent
    binary = ensure_bridge_binary(root)
    threads = str(max(1, min(8, os.cpu_count() or 1)))
    cmd = [str(binary), "-n", "200000000000", "-t", threads, "--no-validate"]
    out = subprocess.check_output(cmd, text=True)
    lines = [line.strip() for line in out.splitlines() if line.strip()]
    if not lines:
        raise RuntimeError("Euler245 bridge produced empty output.")
    return lines[-1]


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

Java

import java.util.*;
import java.math.*;

public class Euler245 {
    static boolean[] isPrimeSmall;
    static int[] primes;

    public static String solve() {
        long LIMIT = 200_000_000_000L;
        int kMax = (int) Math.sqrt((double) LIMIT);
        sieve(kMax);

        long[] rem = new long[kMax + 1];
        int[][] factors = new int[kMax + 1][];
        int[][] exponents = new int[kMax + 1][];
        int[] fcount = new int[kMax + 1];

        List<int[]> tempF = new ArrayList<>();
        List<int[]> tempE = new ArrayList<>();
        for (int k = 2; k <= kMax; k++) {
            rem[k] = (long) k * k - k + 1;
            tempF.add(new int[0]);
            tempE.add(new int[0]);
        }

        for (int p : primes) {
            if (p > kMax)
                break;
            if (p == 2)
                continue;
            if (p == 3) {
                for (int k = 2; k <= kMax; k += 3) {
                    if (rem[k] % 3 != 0)
                        continue;
                    int e = 0;
                    while (rem[k] % 3 == 0) {
                        rem[k] /= 3;
                        e++;
                    }
                    addFactor(k, 3, e, factors, exponents, fcount);
                }
                continue;
            }
            if (p % 3 != 1)
                continue;
            long root = tonelliShanks(p - 3, p);
            long inv2 = (p + 1) / 2;
            long r1 = ((1 + root) % p * inv2) % p;
            long r2 = ((1 + p - root) % p * inv2) % p;
            for (long r : new long[] { r1, r2 }) {
                long kk = r < 2 ? r + ((2 - r + p - 1) / p) * p : r;
                while (kk <= kMax) {
                    int ki = (int) kk;
                    if (rem[ki] % p == 0) {
                        int e = 0;
                        while (rem[ki] % p == 0) {
                            rem[ki] /= p;
                            e++;
                        }
                        addFactor(ki, p, e, factors, exponents, fcount);
                    }
                    kk += p;
                }
            }
        }
        for (int k = 2; k <= kMax; k++) {
            if (rem[k] > 1)
                addFactor(k, (int) rem[k], 1, factors, exponents, fcount);
        }

        Set<Long> results = new HashSet<>();
        for (int k = 2; k <= kMax; k++) {
            long K = (long) k * k - k + 1;
            List<Long> divs = getDivisors(k, factors, exponents, fcount);
            for (long a : divs) {
                long b = K / a;
                if (a > b)
                    continue;
                long pp = k + a, qq = k + b;
                if (pp == qq)
                    continue;
                if (pp > LIMIT / qq)
                    continue;
                if (isPrime(pp) && isPrime(qq))
                    results.add(pp * qq);
            }
        }

        int kRoot = (int) Math.cbrt((double) LIMIT);
        for (int k = 2; k < kRoot; k++) {
            int idx = Arrays.binarySearch(primes, k + 1);
            if (idx < 0)
                idx = -idx - 1;
            dfsMulti(k, idx, 1, 1, 1, 0, LIMIT, results);
        }

        long sum = 0;
        for (long v : results)
            sum += v;
        return String.valueOf(sum);
    }

    static void dfsMulti(int k, int idx, long m, long A, int lastP, int depth, long LIMIT, Set<Long> results) {
        long denom = k * A - ((long) (k - 1)) * m;
        if (denom <= 0)
            return;
        long numer = k * A + 1;
        if (numer % denom == 0) {
            long p = numer / denom;
            if (p > lastP && m <= LIMIT / p && depth >= 2 && isPrime(p))
                results.add(m * p);
        }
        for (int i = idx; i < primes.length; i++) {
            long q = primes[i];
            if (q <= k)
                continue;
            if (m > LIMIT / q / q)
                break;
            dfsMulti(k, i + 1, m * q, A * (q - 1), (int) q, depth + 1, LIMIT, results);
        }
    }

    static int[][] factorsBuf;
    static int[][] exponentsBuf;

    static void addFactor(int k, int p, int e, int[][] factors, int[][] exponents, int[] fcount) {
        int c = fcount[k]++;
        if (factors[k] == null || c >= factors[k].length) {
            int nl = Math.max(4, c + 1);
            factors[k] = Arrays.copyOf(factors[k] == null ? new int[0] : factors[k], nl);
            exponents[k] = Arrays.copyOf(exponents[k] == null ? new int[0] : exponents[k], nl);
        }
        factors[k][c] = p;
        exponents[k][c] = e;
    }

    static List<Long> getDivisors(int k, int[][] factors, int[][] exponents, int[] fcount) {
        List<Long> divs = new ArrayList<>();
        divs.add(1L);
        for (int i = 0; i < fcount[k]; i++) {
            int p = factors[k][i], e = exponents[k][i];
            int sz = divs.size();
            long mult = 1;
            for (int j = 0; j < e; j++) {
                mult *= p;
                for (int m = 0; m < sz; m++)
                    divs.add(divs.get(m) * mult);
            }
        }
        return divs;
    }

    static void sieve(int n) {
        isPrimeSmall = new boolean[n + 1];
        Arrays.fill(isPrimeSmall, true);
        isPrimeSmall[0] = isPrimeSmall[1] = false;
        for (int i = 2; (long) i * i <= n; i++)
            if (isPrimeSmall[i])
                for (int j = i * i; j <= n; j += i)
                    isPrimeSmall[j] = false;
        List<Integer> pl = new ArrayList<>();
        for (int i = 2; i <= n; i++)
            if (isPrimeSmall[i])
                pl.add(i);
        primes = pl.stream().mapToInt(Integer::intValue).toArray();
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        if (n < isPrimeSmall.length)
            return isPrimeSmall[(int) n];
        return millerRabin(n);
    }

    static boolean millerRabin(long n) {
        if (n < 2)
            return false;
        for (long p : new long[] { 2, 3, 5, 7, 11, 13 }) {
            if (n % p == 0)
                return n == p;
        }
        long d = n - 1;
        int s = 0;
        while (d % 2 == 0) {
            d /= 2;
            s++;
        }
        for (long a : new long[] { 2, 3, 5, 7, 11, 13 }) {
            if (a % n == 0)
                continue;
            long x = modPow(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean witness = true;
            for (int i = 0; i < s - 1; i++) {
                x = mulMod(x, x, n);
                if (x == n - 1) {
                    witness = false;
                    break;
                }
            }
            if (witness)
                return false;
        }
        return true;
    }

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

    static long modPow(long base, long exp, long mod) {
        return BigInteger.valueOf(base).modPow(BigInteger.valueOf(exp), BigInteger.valueOf(mod)).longValue();
    }

    static long tonelliShanks(long n, long p) {
        if (n == 0)
            return 0;
        if (p == 2)
            return 1;
        if (p % 4 == 3)
            return modPow(n, (p + 1) / 4, p);
        long q = p - 1;
        int s = 0;
        while (q % 2 == 0) {
            q /= 2;
            s++;
        }
        long z = 2;
        while (modPow(z, (p - 1) / 2, p) != p - 1)
            z++;
        long c = modPow(z, q, p), x = modPow(n, (q + 1) / 2, p), t = modPow(n, q, p);
        int m = s;
        while (t != 1) {
            int i = 1;
            long t2i = mulMod(t, t, p);
            while (t2i != 1) {
                t2i = mulMod(t2i, t2i, p);
                i++;
            }
            long b = modPow(c, 1L << (m - i - 1), p);
            x = mulMod(x, b, p);
            long b2 = mulMod(b, b, p);
            t = mulMod(t, b2, p);
            c = b2;
            m = i;
        }
        return x;
    }

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