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

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$. Every valid composite must be squarefree:

$$n=p_1p_2\cdots p_r,\qquad \varphi(n)=\prod_{i=1}^{r}(p_i-1).$$

This is why the multi-prime branch of the implementation only builds products of distinct primes.

Every prime factor must be larger than $k$

From $k\varphi(n)=(k-1)n+1$ we obtain

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

Because $n$ is squarefree,

$$\frac{n}{\varphi(n)}=\prod_{p\mid n}\frac{p}{p-1}.$$

The function $x\mapsto \frac{x}{x-1}$ decreases for $x\gt 1$. If some prime divisor satisfied $p\le k$, then

$$\frac{p}{p-1}\ge \frac{k}{k-1},$$

and multiplying by the remaining factors, all of which are $\gt 1$, would force $\frac{n}{\varphi(n)}$ to be at least $\frac{k}{k-1}$, contradicting the previous inequality. Hence

$$p_i\gt k.$$

This immediately explains the two main search bounds. If $n=pq$, then $n\gt (k+1)^2$, so only $k\le \sqrt{L}$ can contribute. If $n$ has at least three prime factors, then $n\gt (k+1)^3$, so only $k\lt \sqrt[3]{L}$ need be explored in that branch.

The semiprime case becomes a divisor-pair problem

Let $n=pq$ with distinct primes $p$ and $q$. Then

$$\varphi(n)=(p-1)(q-1)=pq-p-q+1,$$

$$n-\varphi(n)=pq-(pq-p-q+1)=p+q-1.$$

Substituting into $n-1=k(n-\varphi(n))$ gives

$$pq-1=k(p+q-1).$$

Rearranging yields

$$pq-kp-kq+k^2=k^2-k+1,$$

hence

$$\boxed{(p-k)(q-k)=k^2-k+1.}$$

Now define

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

Every semiprime solution for that $k$ comes from a divisor pair $ab=K(k)$:

$$p=k+a,\qquad q=k+b.$$

So the semiprime stage reduces to factoring $K(k)$, enumerating its divisor pairs, shifting them by $k$, and testing primality.

Sieving the polynomial $k^2-k+1$

Factoring each $K(k)$ independently would still be too expensive. The implementations instead sieve the polynomial values arithmetically. If an odd prime $r$ divides $K(k)$, then

$$k^2-k+1\equiv 0\pmod r\iff (2k-1)^2\equiv -3\pmod r.$$

Therefore $k$ must lie in one of the residue classes

$$k\equiv \frac{1\pm \sqrt{-3}}{2}\pmod r.$$

For $r=3$, this collapses to the single class $k\equiv 2\pmod 3$. For other odd primes, a square root of $-3$ exists only when the prime is in the right quadratic-residue class, which here means $r\equiv 1\pmod 3$. That is why the sieve skips $2$ and all primes $r\equiv 2\pmod 3$, and uses Tonelli-Shanks only for the relevant primes to recover the admissible residue classes. Once those classes are known, the prime powers are stripped from every matching $K(k)$, producing the factor tables used later to generate all divisor pairs.

For three or more primes, the last prime is determined explicitly

Now consider a squarefree candidate with at least three prime factors. Suppose we have already chosen distinct primes $q_1,\dots,q_t$ and set

$$m=\prod_{i=1}^{t}q_i,\qquad A=\prod_{i=1}^{t}(q_i-1).$$

If the final number is $n=mp$, then $\varphi(n)=A(p-1)$. Substituting this into

$$k\varphi(n)=(k-1)n+1$$

gives

$$kA(p-1)=(k-1)mp+1.$$

Collecting the terms containing $p$,

$$p\bigl(kA-(k-1)m\bigr)=kA+1,$$

so the final prime is forced to be

$$\boxed{p=\frac{kA+1}{kA-(k-1)m}.}$$

This formula is the backbone of the recursive search. A branch succeeds only if the denominator is positive, it divides the numerator exactly, the resulting $p$ is prime, it is larger than the primes already chosen, and $mp\le L$.

Worked examples

The semiprime $n=15=3\cdot 5$ comes from $k=2$. Here

$$K(2)=2^2-2+1=3,$$

and the divisor pair $(1,3)$ yields

$$p=2+1=3,\qquad q=2+3=5.$$

Indeed $\varphi(15)=8$, so

$$\frac{15-1}{15-\varphi(15)}=\frac{14}{7}=2.$$

A three-prime example is $n=255=3\cdot 5\cdot 17$, again with $k=2$. After choosing $m=3\cdot 5=15$ and $A=(3-1)(5-1)=8$, the forced last prime is

$$p=\frac{2\cdot 8+1}{2\cdot 8-(2-1)\cdot 15}=\frac{17}{1}=17.$$

Since $\varphi(255)=128$, we get

$$\frac{255-1}{255-\varphi(255)}=\frac{254}{127}=2.$$

These two examples are exactly the two mechanisms used by the code: divisor pairs for semiprimes, and the forced-final-prime formula for the squarefree recursive branch.

How the Code Works

Prime tables and primality checks

The full implementations first generate all primes up to $\sqrt{L}$. Those primes drive the congruence sieve for $K(k)$, provide the candidates used in the recursive squarefree search, and answer small primality queries directly. Larger candidates are checked with a deterministic Miller-Rabin test in the relevant 64-bit range.

Semiprime enumeration

For each $k\le \sqrt{L}$, the implementation stores the current remainder of $K(k)$ together with the prime exponents found so far. Each admissible sieve prime contributes one or two residue classes of $k$, and the algorithm walks exactly those progressions while stripping powers of that prime from the corresponding $K(k)$. After the sieve phase, every divisor of $K(k)$ is generated from the recorded factorization, translated into $p=k+a$ and $q=k+b$, tested for primality, and kept only when $pq\le L$.

Recursive squarefree search

The multi-prime branch runs only for $k\lt \sqrt[3]{L}$. It performs a depth-first search over strictly increasing primes $\gt k$, maintaining the pair $(m,A)$. At every node, the formula

$$p=\frac{kA+1}{kA-(k-1)m}$$

is evaluated as a potential closing step. Whenever it produces a valid prime and at least two primes have already been chosen, a composite with at least three factors is recorded. The pruning is strong: if $kA-(k-1)m\le 0$, no completion is possible, and if the next prime $q$ would already force $mq^2\gt L$, then there is no room for both $q$ and a larger final prime.

What the three language implementations do

The C++ implementation contains the full optimized search and parallelizes both collection stages with worker threads. It also includes small-limit brute-force checks to validate the fast method. The Java implementation mirrors the same number-theoretic search directly in one process. The Python implementation is deliberately thin: it ensures a compiled solver is available and then invokes that same fast search, so all three language entries correspond to the same mathematics and the same final answer.

Complexity Analysis

The semiprime stage is the dominant part of the runtime. It covers all $k\le \sqrt{L}$, but its cost is much lower than naive independent factorization because each sieve prime only visits the residue classes where $k^2-k+1$ can actually be divisible by that prime. Divisor enumeration is then done from compact factorizations rather than from trial division.

The multi-prime branch is smaller. It only runs up to $k\lt \sqrt[3]{L}$, only explores increasing primes $\gt k$, and is pruned by the positivity test for $kA-(k-1)m$, by exact divisibility of the forced-final-prime formula, by primality of the recovered last prime, and by the product bound $n\le L$. Memory usage is moderate: a prime list, factor tables for the polynomial values, and a set of accepted solutions before sorting and deduplication. In practice the method is fast because it replaces an impossible totient sweep up to $2\cdot 10^{11}$ with two highly structured searches.

Footnotes and References

Problem page: https://projecteuler.net/problem=245
Euler totient function: Wikipedia - Euler's totient function
Squarefree integer: Wikipedia - Squarefree integer
Tonelli-Shanks algorithm: Wikipedia - Tonelli-Shanks algorithm
Miller-Rabin primality test: Wikipedia - Miller-Rabin primality test

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());
    }
}