Problem 415: Titanic Sets

Project Euler Problem 415 Solution

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

Problem Summary Let $G_N=\{(x,y)\in\mathbb{Z}^2:0\le x,y\le N\}$. A finite set $S\subseteq G_N$ is called titanic if some line passes through exactly two points of $S$. We must compute $T(N)$, the number of titanic subsets of $G_N$, modulo $10^8$. The implementation never enumerates subsets directly. Instead, it counts the complement: empty sets, singletons, and larger subsets whose selected points are all collinear. Mathematical Approach Step 1: Reduce to the Complement The key theorem is Sylvester-Gallai: every finite non-collinear set of real points has an ordinary line, meaning a line containing exactly two of the points. Therefore a lattice subset is not titanic if and only if it is empty, a singleton, or all of its points lie on one common line. If $\operatorname{NT}(N)$ denotes the number of non-titanic subsets, then $$T(N)=2^{(N+1)^2}-\operatorname{NT}(N)\pmod{10^8}.$$ Step 2: Replace Subset Counting by Line-Length Counting Suppose a grid line contains exactly $k$ lattice points. The non-titanic subsets supported on that line are exactly the subsets of size at least $3$, so that line contributes $$\sum_{r=3}^{k}\binom{k}{r}=2^k-1-k-\binom{k}{2}.$$ Now define $F_\ell(N)$ as the number of grid lines containing at least $\ell+1$ lattice points....

Detailed mathematical approach

Problem Summary

Let $G_N=\{(x,y)\in\mathbb{Z}^2:0\le x,y\le N\}$. A finite set $S\subseteq G_N$ is called titanic if some line passes through exactly two points of $S$. We must compute $T(N)$, the number of titanic subsets of $G_N$, modulo $10^8$.

The implementation never enumerates subsets directly. Instead, it counts the complement: empty sets, singletons, and larger subsets whose selected points are all collinear.

Mathematical Approach

Step 1: Reduce to the Complement

The key theorem is Sylvester-Gallai: every finite non-collinear set of real points has an ordinary line, meaning a line containing exactly two of the points. Therefore a lattice subset is not titanic if and only if it is empty, a singleton, or all of its points lie on one common line.

If $\operatorname{NT}(N)$ denotes the number of non-titanic subsets, then

$$T(N)=2^{(N+1)^2}-\operatorname{NT}(N)\pmod{10^8}.$$

Step 2: Replace Subset Counting by Line-Length Counting

Suppose a grid line contains exactly $k$ lattice points. The non-titanic subsets supported on that line are exactly the subsets of size at least $3$, so that line contributes

$$\sum_{r=3}^{k}\binom{k}{r}=2^k-1-k-\binom{k}{2}.$$

Now define $F_\ell(N)$ as the number of grid lines containing at least $\ell+1$ lattice points. For a fixed line with $k$ points, one has the identity

$$\sum_{\ell=2}^{k-1}\bigl(2^\ell-(\ell+1)\bigr)=2^k-1-k-\binom{k}{2}.$$

Therefore

$$\operatorname{NT}(N)=1+(N+1)^2+\sum_{\ell=2}^{N}F_\ell(N)\bigl(2^\ell-(\ell+1)\bigr).$$

The whole task is now reduced to computing $F_\ell(N)$ efficiently for all $\ell\ge2$.

Step 3: Primitive Directions and Segment Placements

Write $a=N+1$. Every non-axis grid line has a primitive step vector $(u,v)$ with $1\le v\le u$ and $\gcd(u,v)=1$. The pair $(u,v)$ represents the slope classes $\pm v/u$, and when $u\ne v$, also the swapped classes $\pm u/v$. Horizontal and vertical lines are handled separately.

For such a primitive direction, the number of placements of an $\ell$-step segment inside the square is

$$C_\ell(u,v)=\max(a-\ell u,0)\max(a-\ell v,0).$$

If a line contains $k$ lattice points, then it contributes exactly $k-\ell$ segments of length $\ell$ and $k-\ell-1$ segments of length $\ell+1$. Their difference is $1$ exactly when $k\ge\ell+1$. Hence the number of distinct lines in this direction class that contain at least $\ell+1$ points is

$$C_\ell(u,v)-C_{\ell+1}(u,v).$$

After doubling for the two signs and adding the $a$ horizontal plus $a$ vertical lines, we recover $F_\ell(N)$.

Step 4: Totient-Based Direction Statistics

Let

$$\Phi(M)=\sum_{d\le M}\varphi(d),\qquad \Psi_1(M)=\sum_{d\le M} d\,\varphi(d),\qquad \Psi_2(M)=\sum_{d\le M} d^2\varphi(d).$$

The implementation stores the primitive-direction information in three cumulative quantities:

$$c(M)=2\Phi(M)-1,\qquad s_{xy}(M)=3\Psi_1(M)-1,\qquad s_{\prod}(M)=\Psi_2(M).$$

Why do these formulas appear? Group primitive directions by their larger coordinate $d$. For fixed $d>1$, there are $\varphi(d)$ reduced numerators, their sum is $d\varphi(d)/2$, and the swapped directions contribute the same amount. The exceptional diagonal case $d=1$ explains the subtractive constants.

Step 5: Fast Evaluation of $\Phi$, $\Psi_1$, and $\Psi_2$

The code does not recompute these sums from scratch for every query. Using $\sum_{d\mid n}\varphi(d)=n$, one gets the divisor-summatory identities

$$\sum_{q=1}^{n}\Phi\!\left(\left\lfloor\frac{n}{q}\right\rfloor\right)=\sum_{m=1}^{n}m=\frac{n(n+1)}{2},$$

$$\sum_{q=1}^{n}q\,\Psi_1\!\left(\left\lfloor\frac{n}{q}\right\rfloor\right)=\sum_{m=1}^{n}m^2=\frac{n(n+1)(2n+1)}{6},$$

$$\sum_{q=1}^{n}q^2\,\Psi_2\!\left(\left\lfloor\frac{n}{q}\right\rfloor\right)=\sum_{m=1}^{n}m^3=\left(\frac{n(n+1)}{2}\right)^2.$$

Those identities lead to memoized recurrences over quotient blocks $\left\lfloor n/q\right\rfloor$, which is why the summatory functions can be queried many times without a linear scan up to $n$.

Step 6: Closed Form for $F_\ell(N)$

Set

$$M=\left\lfloor\frac{N}{\ell}\right\rfloor,\qquad M'=\left\lfloor\frac{N}{\ell+1}\right\rfloor,\qquad a=N+1.$$

On any interval of $\ell$ where $M$ and $M'$ stay constant, the line count becomes a quadratic polynomial in $\ell$:

$$F_\ell(N)=2a+2\left(q_2\ell^2+q_1\ell+q_0\right),$$

where

$$q_2=s_{\prod}(M)-s_{\prod}(M'),$$

$$q_1=a\bigl(s_{xy}(M')-s_{xy}(M)\bigr)-2s_{\prod}(M'),$$

$$q_0=a^2\bigl(c(M)-c(M')\bigr)+a\,s_{xy}(M')-s_{\prod}(M').$$

This is exactly the quadratic block formula that appears in the implementation.

Step 7: Block Summation

The values of $\left\lfloor N/\ell\right\rfloor$ change only $O(\sqrt N)$ times, so the final summation is performed block by block. For $\ell\in[L_1,L_2]$, the contribution is

$$\sum_{\ell=L_1}^{L_2}F_\ell(N)\bigl(2^\ell-(\ell+1)\bigr).$$

After substituting the quadratic form, the implementation only needs closed forms for

$$\sum 2^\ell,\qquad \sum \ell 2^\ell,\qquad \sum \ell^2 2^\ell,\qquad \sum \ell,\qquad \sum \ell^2,\qquad \sum \ell^3.$$

The exponential identities used are

$$\sum_{j=0}^{n}2^j=2^{n+1}-1,\qquad \sum_{j=0}^{n}j2^j=(n-1)2^{n+1}+2,$$

$$\sum_{j=0}^{n}j^2 2^j=(n^2-2n+3)2^{n+1}-6.$$

Worked Checkpoint: $N=4$

For the $5\times5$ grid, the formula gives

$$F_2(4)=32,\qquad F_3(4)=16,\qquad F_4(4)=12.$$

Hence

$$\operatorname{NT}(4)=1+25+32(1)+16(4)+12(11)=254,$$

and therefore

$$T(4)=2^{25}-254=33554178,$$

which matches the checkpoint verified by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same mathematical pipeline. They precompute Euler totients up to a fixed cutoff with a linear sieve, build prefix tables for $\Phi$, $\Psi_1$, and $\Psi_2$, and memoize larger arguments so that each summatory query is solved only once.

After that, the implementation enumerates the distinct quotient blocks of $\lfloor N/\ell\rfloor$. Inside each block it evaluates the quadratic line-count formula, applies the closed forms for the exponential and polynomial sums, accumulates the complement count, and finally subtracts that value from $2^{(N+1)^2}$ modulo $10^8$.

Complexity Analysis

Let $B$ be the totient precomputation cutoff and let $Q$ be the number of distinct values taken by $\lfloor N/\ell\rfloor$; one has $Q=O(\sqrt N)$. The sieve and prefix tables cost near-linear time and $O(B)$ memory. The outer summation uses only $O(Q)$ quotient blocks, and the cached summatory totient queries are shared across those blocks instead of being recomputed. In practice, the heavy work is the prefix-table construction plus the $O(\sqrt N)$-scale block arithmetic.

References

Problem page: https://projecteuler.net/problem=415
Sylvester-Gallai theorem: Wikipedia — Sylvester-Gallai theorem
Euler's totient function: Wikipedia — Euler's totient function
Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method
Graham, Knuth, Patashnik, Concrete Mathematics, 2nd ed., Addison-Wesley, Chapter 2.

Problem 415 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <unordered_map>
#include <vector>

using namespace std;

namespace {

constexpr uint64_t MOD = 100000000ULL;
constexpr long long DEFAULT_N = 100000000000LL;
constexpr long long PRECOMP_LIMIT = 25000000LL;

struct FastHash {
    size_t operator()(uint64_t x) const noexcept {
        static const uint64_t kMul = 0x9e3779b97f4a7c15ULL;
        x += kMul;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
        x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
        return static_cast<size_t>(x ^ (x >> 31));
    }
};

inline uint64_t mod_norm(int64_t v) {
    v %= static_cast<int64_t>(MOD);
    if (v < 0) v += static_cast<int64_t>(MOD);
    return static_cast<uint64_t>(v);
}

inline uint64_t mod_mul(uint64_t a, uint64_t b) {
    return static_cast<uint64_t>((__int128)(a % MOD) * (b % MOD) % MOD);
}

uint64_t pow2_mod(unsigned __int128 exp) {
    uint64_t base = 2 % MOD;
    uint64_t res = 1;
    while (exp > 0) {
        if (exp & 1) res = mod_mul(res, base);
        base = mod_mul(base, base);
        exp >>= 1;
    }
    return res;
}

uint64_t sum1_upto(long long n) {
    if (n <= 0) return 0;
    __int128 v = static_cast<__int128>(n) * (n + 1) / 2;
    return static_cast<uint64_t>(v % MOD);
}

uint64_t sum2_upto(long long n) {
    if (n <= 0) return 0;
    __int128 v = static_cast<__int128>(n) * (n + 1) * (2 * n + 1) / 6;
    return static_cast<uint64_t>(v % MOD);
}

uint64_t sum3_upto(long long n) {
    if (n <= 0) return 0;
    uint64_t a = static_cast<uint64_t>(n % static_cast<long long>(MOD));
    uint64_t b = static_cast<uint64_t>((n + 1) % static_cast<long long>(MOD));
    if ((n & 1LL) == 0) {
        a = static_cast<uint64_t>((n / 2) % static_cast<long long>(MOD));
    } else {
        b = static_cast<uint64_t>(((n + 1) / 2) % static_cast<long long>(MOD));
    }
    uint64_t t = mod_mul(a, b);
    return mod_mul(t, t);
}

uint64_t sum_range_1(long long l, long long r) {
    return mod_norm(static_cast<int64_t>(sum1_upto(r)) - static_cast<int64_t>(sum1_upto(l - 1)));
}

uint64_t sum_range_2(long long l, long long r) {
    return mod_norm(static_cast<int64_t>(sum2_upto(r)) - static_cast<int64_t>(sum2_upto(l - 1)));
}

uint64_t sum_range_3(long long l, long long r) {
    return mod_norm(static_cast<int64_t>(sum3_upto(r)) - static_cast<int64_t>(sum3_upto(l - 1)));
}

uint64_t sum_q_mod(long long l, long long r) {
    __int128 cnt = static_cast<__int128>(r - l + 1);
    __int128 sum = static_cast<__int128>(l + r) * cnt / 2;
    return static_cast<uint64_t>(sum % MOD);
}

uint64_t sum_lpow_upto(long long n, uint64_t pow2_n_plus1) {
    if (n < 0) return 0;
    uint64_t n_mod = static_cast<uint64_t>(n % static_cast<long long>(MOD));
    uint64_t coef = mod_norm(static_cast<int64_t>(n_mod) - 1);
    return (mod_mul(coef, pow2_n_plus1) + 2) % MOD;
}

uint64_t sum_l2pow_upto(long long n, uint64_t pow2_n_plus1) {
    if (n < 0) return 0;
    uint64_t n_mod = static_cast<uint64_t>(n % static_cast<long long>(MOD));
    uint64_t n2_mod = mod_mul(n_mod, n_mod);
    uint64_t coef = mod_norm(static_cast<int64_t>(n2_mod) -
                             static_cast<int64_t>(mod_mul(2, n_mod)) + 3);
    uint64_t val = mod_mul(coef, pow2_n_plus1);
    return mod_norm(static_cast<int64_t>(val) - 6);
}

struct TotientSummatory {
    int limit = 0;
    vector<int> phi;
    vector<uint32_t> pref_phi;
    vector<uint32_t> pref_weighted;
    vector<uint32_t> pref_weighted2;
    unordered_map<long long, uint64_t, FastHash> cache_phi;
    unordered_map<long long, uint64_t, FastHash> cache_weighted;
    unordered_map<long long, uint64_t, FastHash> cache_weighted2;

    explicit TotientSummatory(long long max_n) {
        limit = static_cast<int>(min(max_n, PRECOMP_LIMIT));
        if (limit < 1) limit = 1;
        phi.assign(limit + 1, 0);
        pref_phi.assign(limit + 1, 0);
        pref_weighted.assign(limit + 1, 0);
        pref_weighted2.assign(limit + 1, 0);

        vector<int> primes;
        vector<uint8_t> is_comp(limit + 1, 0);
        phi[1] = 1;
        for (int i = 2; i <= limit; ++i) {
            if (!is_comp[i]) {
                primes.push_back(i);
                phi[i] = i - 1;
            }
            for (int p : primes) {
                int64_t v = 1LL * i * p;
                if (v > limit) break;
                is_comp[static_cast<int>(v)] = 1;
                if (i % p == 0) {
                    phi[static_cast<int>(v)] = phi[i] * p;
                    break;
                }
                phi[static_cast<int>(v)] = phi[i] * (p - 1);
            }
        }

        for (int i = 1; i <= limit; ++i) {
            uint64_t add_phi = static_cast<uint64_t>(phi[i]) % MOD;
            pref_phi[i] = static_cast<uint32_t>((pref_phi[i - 1] + add_phi) % MOD);

            uint64_t add_weight = (static_cast<uint64_t>(i) * phi[i]) % MOD;
            pref_weighted[i] = static_cast<uint32_t>((pref_weighted[i - 1] + add_weight) % MOD);

            uint64_t i_mod = static_cast<uint64_t>(i) % MOD;
            uint64_t i2_mod = mod_mul(i_mod, i_mod);
            uint64_t add_weight2 = mod_mul(i2_mod, static_cast<uint64_t>(phi[i]) % MOD);
            pref_weighted2[i] = static_cast<uint32_t>((pref_weighted2[i - 1] + add_weight2) % MOD);
        }

        cache_phi.reserve(1 << 20);
        cache_weighted.reserve(1 << 20);
        cache_weighted2.reserve(1 << 20);
    }

    uint64_t phi_sum(long long n) {
        if (n <= limit) return pref_phi[static_cast<int>(n)];
        auto it = cache_phi.find(n);
        if (it != cache_phi.end()) return it->second;

        uint64_t res = sum1_upto(n);
        for (long long l = 2; l <= n; ) {
            long long t = n / l;
            long long r = n / t;
            uint64_t sub = mod_mul(static_cast<uint64_t>(r - l + 1), phi_sum(t));
            res = (res + MOD - sub) % MOD;
            l = r + 1;
        }
        cache_phi[n] = res;
        return res;
    }

    uint64_t phi_weighted_sum(long long n) {
        if (n <= limit) return pref_weighted[static_cast<int>(n)];
        auto it = cache_weighted.find(n);
        if (it != cache_weighted.end()) return it->second;

        uint64_t res = sum2_upto(n);
        for (long long l = 2; l <= n; ) {
            long long t = n / l;
            long long r = n / t;
            uint64_t sum_q = sum_q_mod(l, r);
            uint64_t sub = mod_mul(sum_q, phi_weighted_sum(t));
            res = (res + MOD - sub) % MOD;
            l = r + 1;
        }
        cache_weighted[n] = res;
        return res;
    }

    uint64_t phi_weighted2_sum(long long n) {
        if (n <= limit) return pref_weighted2[static_cast<int>(n)];
        auto it = cache_weighted2.find(n);
        if (it != cache_weighted2.end()) return it->second;

        uint64_t res = sum3_upto(n);
        for (long long l = 2; l <= n; ) {
            long long t = n / l;
            long long r = n / t;
            uint64_t sum_q2 = sum_range_2(l, r);
            uint64_t sub = mod_mul(sum_q2, phi_weighted2_sum(t));
            res = (res + MOD - sub) % MOD;
            l = r + 1;
        }
        cache_weighted2[n] = res;
        return res;
    }
};

struct RawRange {
    long long L1;
    long long L2;
    long long M;
    long long M2;
};

struct Range {
    long long L1;
    long long L2;
    uint64_t c1;
    uint64_t sxy1;
    uint64_t sprod1;
    uint64_t c2;
    uint64_t sxy2;
    uint64_t sprod2;
};

uint64_t compute_T(long long N, TotientSummatory& ts, unsigned threads) {
    unsigned __int128 exp = static_cast<unsigned __int128>(N + 1) * (N + 1);
    uint64_t total_subsets = pow2_mod(exp);

    uint64_t points_mod = mod_mul(static_cast<uint64_t>((N + 1) % MOD),
                                  static_cast<uint64_t>((N + 1) % MOD));

    if (N < 2) {
        uint64_t non_titanic = (1 + points_mod) % MOD;
        return mod_norm(static_cast<int64_t>(total_subsets) - static_cast<int64_t>(non_titanic));
    }

    vector<RawRange> raw_ranges;
    vector<long long> ms;
    for (long long L = 2; L <= N; ) {
        long long M1 = N / L;
        long long M2 = (L + 1 <= N) ? (N / (L + 1)) : 0;
        long long R1 = N / M1;
        long long R2 = (M2 == 0) ? N : (N / M2 - 1);
        long long R = min(R1, R2);
        raw_ranges.push_back({L, R, M1, M2});
        ms.push_back(M1);
        if (M2 > 0) ms.push_back(M2);
        L = R + 1;
    }

    sort(ms.begin(), ms.end());
    ms.erase(unique(ms.begin(), ms.end()), ms.end());

    struct PrecompVals {
        uint64_t c;
        uint64_t sxy;
        uint64_t sprod;
    };

    unordered_map<long long, PrecompVals, FastHash> precomp;
    precomp.reserve(ms.size() * 2);
    for (long long M : ms) {
        uint64_t sum_phi = ts.phi_sum(M);
        uint64_t sum_weighted = ts.phi_weighted_sum(M);
        uint64_t sum_weighted2 = ts.phi_weighted2_sum(M);
        uint64_t c = (2 * sum_phi + MOD - 1) % MOD;
        uint64_t sxy = (3 * sum_weighted + MOD - 1) % MOD;
        uint64_t sprod = sum_weighted2 % MOD;
        precomp.emplace(M, PrecompVals{c, sxy, sprod});
    }

    vector<Range> ranges;
    ranges.reserve(raw_ranges.size());
    for (const auto& rr : raw_ranges) {
        auto it1 = precomp.find(rr.M);
        PrecompVals v1 = it1->second;
        PrecompVals v2{0, 0, 0};
        if (rr.M2 > 0) {
            auto it2 = precomp.find(rr.M2);
            v2 = it2->second;
        }
        ranges.push_back({rr.L1, rr.L2, v1.c, v1.sxy, v1.sprod, v2.c, v2.sxy, v2.sprod});
    }

    uint64_t a = static_cast<uint64_t>((N + 1) % static_cast<long long>(MOD));
    uint64_t a2 = mod_mul(a, a);

    if (threads == 0) threads = 1;
    if (ranges.size() < 2000) threads = 1;
    if (threads > ranges.size()) threads = static_cast<unsigned>(ranges.size());

    vector<uint64_t> partial(threads, 0);
    atomic<size_t> next{0};
    const size_t chunk = 128;

    auto worker = [&](unsigned idx) {
        uint64_t local = 0;
        while (true) {
            size_t start = next.fetch_add(chunk, memory_order_relaxed);
            if (start >= ranges.size()) break;
            size_t end = min(start + chunk, ranges.size());
            for (size_t i = start; i < end; ++i) {
                const Range& rg = ranges[i];
                uint64_t q2 = mod_norm(static_cast<int64_t>(rg.sprod1) - static_cast<int64_t>(rg.sprod2));
                uint64_t q1 = mod_norm(static_cast<int64_t>(mod_norm(
                                     static_cast<int64_t>(mod_mul(a, rg.sxy2)) -
                                     static_cast<int64_t>(mod_mul(a, rg.sxy1))) -
                                     static_cast<int64_t>(mod_mul(2, rg.sprod2))));
                uint64_t q0 = mod_norm(static_cast<int64_t>(mod_mul(a2,
                                     mod_norm(static_cast<int64_t>(rg.c1) - static_cast<int64_t>(rg.c2)))) +
                                     static_cast<int64_t>(mod_mul(a, rg.sxy2)) -
                                     static_cast<int64_t>(rg.sprod2));

                uint64_t p2 = mod_mul(2, q2);
                uint64_t p1 = mod_mul(2, q1);
                uint64_t p0 = mod_norm(static_cast<int64_t>(mod_mul(2, q0)) +
                                       static_cast<int64_t>(mod_mul(2, a)));

                long long L1 = rg.L1;
                long long L2 = rg.L2;
                uint64_t pow_L1 = pow2_mod(static_cast<unsigned __int128>(L1));
                uint64_t pow_L2p1 = pow2_mod(static_cast<unsigned __int128>(L2) + 1);
                uint64_t s0_r = mod_norm(static_cast<int64_t>(pow_L2p1) - 1);
                uint64_t s0_l = mod_norm(static_cast<int64_t>(pow_L1) - 1);
                uint64_t s0 = mod_norm(static_cast<int64_t>(s0_r) - static_cast<int64_t>(s0_l));

                uint64_t s1_r = sum_lpow_upto(L2, pow_L2p1);
                uint64_t s1_l = sum_lpow_upto(L1 - 1, pow_L1);
                uint64_t s1 = mod_norm(static_cast<int64_t>(s1_r) - static_cast<int64_t>(s1_l));

                uint64_t s2_r = sum_l2pow_upto(L2, pow_L2p1);
                uint64_t s2_l = sum_l2pow_upto(L1 - 1, pow_L1);
                uint64_t s2 = mod_norm(static_cast<int64_t>(s2_r) - static_cast<int64_t>(s2_l));

                uint64_t t1 = sum_range_1(L1, L2);
                uint64_t t2 = sum_range_2(L1, L2);
                uint64_t t3 = sum_range_3(L1, L2);
                uint64_t cnt = static_cast<uint64_t>((L2 - L1 + 1) % MOD);

                uint64_t sum_pow = 0;
                sum_pow = (sum_pow + mod_mul(p2, s2)) % MOD;
                sum_pow = (sum_pow + mod_mul(p1, s1)) % MOD;
                sum_pow = (sum_pow + mod_mul(p0, s0)) % MOD;

                uint64_t sum_poly = 0;
                uint64_t t3_plus_t2 = (t3 + t2) % MOD;
                uint64_t t2_plus_t1 = (t2 + t1) % MOD;
                uint64_t t1_plus_cnt = (t1 + cnt) % MOD;
                sum_poly = (sum_poly + mod_mul(p2, t3_plus_t2)) % MOD;
                sum_poly = (sum_poly + mod_mul(p1, t2_plus_t1)) % MOD;
                sum_poly = (sum_poly + mod_mul(p0, t1_plus_cnt)) % MOD;

                uint64_t term = mod_norm(static_cast<int64_t>(sum_pow) - static_cast<int64_t>(sum_poly));

                local += term;
                if (local >= MOD) local %= MOD;
            }
        }
        partial[idx] = local % MOD;
    };

    if (threads == 1) {
        worker(0);
    } else {
        vector<thread> pool;
        pool.reserve(threads);
        for (unsigned t = 0; t < threads; ++t) pool.emplace_back(worker, t);
        for (auto& th : pool) th.join();
    }

    uint64_t sum_line = 0;
    for (uint64_t v : partial) {
        sum_line += v;
        if (sum_line >= MOD) sum_line %= MOD;
    }

    uint64_t non_titanic = (1 + points_mod + sum_line) % MOD;
    return mod_norm(static_cast<int64_t>(total_subsets) - static_cast<int64_t>(non_titanic));
}

bool validate(TotientSummatory& ts) {
    struct Test {
        long long n;
        uint64_t expected;
    };
    const Test tests[] = {
        {1, 11},
        {2, 494},
        {4, 33554178},
        {111, 13500401},
        {100000, 63259062},
    };

    for (const auto& test : tests) {
        uint64_t got = compute_T(test.n, ts, 1);
        if (got != test.expected) {
            cerr << "Validation failed for N=" << test.n
                 << ": got " << got << ", expected " << test.expected << '\n';
            return false;
        }
    }
    return true;
}

} // namespace

int main(int argc, char** argv) {
    long long target = DEFAULT_N;
    if (argc > 1) {
        target = max(0LL, atoll(argv[1]));
    }
    unsigned threads = thread::hardware_concurrency();
    if (argc > 2) {
        threads = max(1, atoi(argv[2]));
    }

    long long max_n = max(target, 100000LL);
    TotientSummatory ts(max_n);
    if (!validate(ts)) return 1;

    uint64_t result = compute_T(target, ts, threads);
    cout << result << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler415 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler415.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler415_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler415 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler415 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler415 C++ bridge produced empty output.");
        }
        return parsed;
    }

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