Problem 735: Divisors of $2n^2$

Project Euler Problem 735 Solution

EulerSolve provides an optimized solution for Project Euler Problem 735, Divisors of $2n^2$, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $n$, define $$a(n)=\#\left\{d\in\mathbb{Z}_{>0}: d\mid 2n^2,\ d\le n\right\}.$$ The problem asks for $$F(N)=\sum_{n=1}^{N} a(n).$$ A brute-force method would test every $d\le n$ for every $n\le N$, which is far too slow for the large target value. The successful approach is to count admissible pairs $(n,d)$ in a structured way, then compress the remaining sums with Möbius inversion and quotient-block floor sums. Mathematical Approach Instead of viewing the problem as a separate divisor search for each $n$, we count all valid pairs $(n,d)$ at once. Step 1: Recast the Problem as a Pair Count By definition, $F(N)$ counts pairs $(n,d)$ such that $$1\le n\le N,\qquad d\mid 2n^2,\qquad d\le n.$$ So we may write $$F(N)=\#\left\{(n,d):1\le n\le N,\ d\mid 2n^2,\ d\le n\right\}.$$ The inequality $d\le n$ and the divisibility condition interact naturally with $\gcd(n,d)$, so that is the right quantity to isolate first....

Detailed mathematical approach

Problem Summary

For each positive integer $n$, define

$$a(n)=\#\left\{d\in\mathbb{Z}_{>0}: d\mid 2n^2,\ d\le n\right\}.$$

The problem asks for

$$F(N)=\sum_{n=1}^{N} a(n).$$

A brute-force method would test every $d\le n$ for every $n\le N$, which is far too slow for the large target value. The successful approach is to count admissible pairs $(n,d)$ in a structured way, then compress the remaining sums with Möbius inversion and quotient-block floor sums.

Mathematical Approach

Instead of viewing the problem as a separate divisor search for each $n$, we count all valid pairs $(n,d)$ at once.

Step 1: Recast the Problem as a Pair Count

By definition, $F(N)$ counts pairs $(n,d)$ such that

$$1\le n\le N,\qquad d\mid 2n^2,\qquad d\le n.$$

So we may write

$$F(N)=\#\left\{(n,d):1\le n\le N,\ d\mid 2n^2,\ d\le n\right\}.$$

The inequality $d\le n$ and the divisibility condition interact naturally with $\gcd(n,d)$, so that is the right quantity to isolate first.

Step 2: Split Off the Common GCD

For any counted pair, set

$$g=\gcd(n,d),\qquad n=g\,m,\qquad d=g\,u,\qquad \gcd(m,u)=1.$$

Because $d\le n$, we immediately have

$$u\le m.$$

Also, from $d\mid 2n^2$ we get

$$g\,u \mid 2g^2m^2,$$

and after dividing by $g$,

$$u\mid 2g\,m^2.$$

Since $\gcd(u,m)=1$, every prime factor of $u$ must come from $2g$, so in fact

$$u\mid 2g.$$

That condition leads to two natural cases: $u$ odd and $u$ even.

Step 3: The Odd Branch

If $u$ is odd, then $u\mid 2g$ forces $u\mid g$. Write

$$g=u\,t.$$

Substituting back gives

$$d=u^2t,\qquad n=u\,m\,t,\qquad \gcd(m,u)=1,\qquad m\ge u.$$

For fixed $u$ and $m$, the parameter $t$ can be any positive integer with $u m t\le N$, so it contributes

$$\left\lfloor\frac{N}{u m}\right\rfloor$$

choices. Therefore the odd contribution is

$$C_{\mathrm{odd}}(N)=\sum_{\substack{u\le \sqrt N\\u\text{ odd}}}\ \sum_{\substack{m\ge u\\ \gcd(m,u)=1}}\left\lfloor\frac{N}{u m}\right\rfloor.$$

The range $u\le \sqrt N$ is automatic because $n=u m t\ge u^2$.

Step 4: The Even Branch

If $u$ is even, write

$$u=2w.$$

Then $u\mid 2g$ becomes $w\mid g$, so we may set

$$g=w\,t.$$

Now

$$d=2w^2t,\qquad n=w\,m\,t,\qquad \gcd(m,2w)=1,\qquad m\ge 2w.$$

The condition $\gcd(m,2w)=1$ means in particular that $m$ must be odd, and also $\gcd(m,w)=1$. For fixed $w$ and $m$, the number of possible $t$ is

$$\left\lfloor\frac{N}{w m}\right\rfloor.$$

So the even contribution is

$$C_{\mathrm{even}}(N)=\sum_{w\le \sqrt{N/2}}\ \sum_{\substack{m\ge 2w\\ m\text{ odd}\\ \gcd(m,w)=1}}\left\lfloor\frac{N}{w m}\right\rfloor.$$

Again the range follows from $n=wmt\ge 2w^2$. The full answer is simply

$$F(N)=C_{\mathrm{odd}}(N)+C_{\mathrm{even}}(N).$$

Step 5: Remove Coprimality with Möbius Inversion

The remaining obstacle is the condition $\gcd(m,u)=1$ or $\gcd(m,w)=1$. This is handled with the Möbius function:

$$[\gcd(x,y)=1]=\sum_{d\mid \gcd(x,y)} \mu(d).$$

Applying this to the odd branch gives

$$C_{\mathrm{odd}}(N)=\sum_{\substack{u\le \sqrt N\\u\text{ odd}}}\ \sum_{d\mid u}^{\mathrm{sqfree}}\mu(d)\sum_{r\ge \lceil u/d\rceil}\left\lfloor\frac{\left\lfloor N/(ud)\right\rfloor}{r}\right\rfloor.$$

Here $d$ runs over the squarefree divisors of $u$, and we rewrote $m=d r$.

For the even branch, only odd squarefree divisors matter because $m$ is already odd. That yields

$$C_{\mathrm{even}}(N)=\sum_{w\le \sqrt{N/2}}\ \sum_{\substack{d\mid w\\ d\text{ odd}}}^{\mathrm{sqfree}}\mu(d)\sum_{\substack{r\ge \lceil 2w/d\rceil\\ r\text{ odd}}}\left\lfloor\frac{\left\lfloor N/(wd)\right\rfloor}{r}\right\rfloor.$$

These are exactly the two families of sums evaluated by the implementation.

Step 6: Convert the Odd-Denominator Sum into Two Standard Floor Sums

Define the tail floor sum

$$\operatorname{FS}(X,L)=\sum_{m=L}^{X}\left\lfloor\frac{X}{m}\right\rfloor.$$

Then the odd-denominator version is obtained by subtracting the even denominators:

$$\sum_{\substack{m\ge L\\ m\text{ odd}}}\left\lfloor\frac{X}{m}\right\rfloor=\operatorname{FS}(X,L)-\operatorname{FS}\left(\left\lfloor\frac{X}{2}\right\rfloor,\left\lceil\frac{L}{2}\right\rceil\right).$$

This identity is the key to turning the even branch into the same kind of floor-sum object as the odd branch. The implementation then evaluates $\operatorname{FS}(X,L)$ in quotient blocks, using the fact that $\left\lfloor X/m\right\rfloor$ is constant on long intervals.

Worked Example: $F(15)=63$

The implementation checks the small value $F(15)=63$. Using the branch formulas, we can see why.

For the odd branch, $u\le \sqrt{15}$, so $u=1$ or $u=3$.

When $u=1$, every $m\ge 1$ is coprime to $1$, so

$$\sum_{m=1}^{15}\left\lfloor\frac{15}{m}\right\rfloor=45.$$

When $u=3$, we need $m\ge 3$, $\gcd(m,3)=1$, and $\left\lfloor 15/(3m)\right\rfloor>0$. Only $m=4,5$ contribute, giving

$$\left\lfloor\frac{15}{12}\right\rfloor+\left\lfloor\frac{15}{15}\right\rfloor=1+1=2.$$

Hence

$$C_{\mathrm{odd}}(15)=45+2=47.$$

For the even branch, $w\le \sqrt{15/2}$, so $w=1,2$.

When $w=1$, the odd integers $m\ge 2$ that contribute are $3,5,7,9,11,13,15$, so

$$\left\lfloor\frac{15}{3}\right\rfloor+\left\lfloor\frac{15}{5}\right\rfloor+\left\lfloor\frac{15}{7}\right\rfloor+\left\lfloor\frac{15}{9}\right\rfloor+\left\lfloor\frac{15}{11}\right\rfloor+\left\lfloor\frac{15}{13}\right\rfloor+\left\lfloor\frac{15}{15}\right\rfloor=14.$$

When $w=2$, we need odd $m\ge 4$, and only $m=5,7$ contribute:

$$\left\lfloor\frac{15}{10}\right\rfloor+\left\lfloor\frac{15}{14}\right\rfloor=1+1=2.$$

Thus

$$C_{\mathrm{even}}(15)=14+2=16,$$

and finally

$$F(15)=47+16=63.$$

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical decomposition. The optimized computation first determines the two outer ranges, $\lfloor\sqrt N\rfloor$ for the odd branch and $\lfloor\sqrt{N/2}\rfloor$ for the even branch.

Next, it precomputes for every integer up to that range all squarefree divisors together with the corresponding Möbius sign. This is done by extracting the distinct prime factors of each integer and enumerating their subsets, because every subset produces one squarefree divisor and its sign $(-1)^k$.

The main summation then iterates through the odd and even branches exactly as in the formulas above. Each inner tail sum is reduced to $\operatorname{FS}(X,L)$, and $\operatorname{FS}(X,L)$ is not evaluated term-by-term. Instead, the implementation groups consecutive denominators that give the same quotient $\left\lfloor X/m\right\rfloor$, so one arithmetic step handles an entire block.

For large inputs, the C++ implementation can split the odd and even ranges into independent blocks and process them in parallel. The Python and Java implementations reuse that same optimized native computation and only return the final numeric output. Before the full run, the program also checks small checkpoints and compares against a direct brute-force count on a small limit.

Complexity Analysis

Let

$$M=\max\left(\left\lfloor\sqrt N\right\rfloor,\left\lfloor\sqrt{N/2}\right\rfloor\right).$$

Building the smallest-prime-factor sieve up to $M$ costs $O(M\log\log M)$ time and $O(M)$ memory. The stored squarefree-divisor cache contains $\sum_{k\le M}2^{\omega(k)}$ entries, whose average order is $O(M\log M)$, so that cache dominates memory usage in practice.

The two outer summations only run to order $\sqrt N$, not to $N$. For each cached divisor entry, the remaining work is a quotient-block floor sum rather than a linear scan over all denominators. The exact worst-case expression is messy because it depends on both divisor structure and quotient-block lengths, but the practical behavior is close to $N^{1/2+o(1)}$, which is dramatically better than the quadratic brute-force definition.

Footnotes and References

Problem page: https://projecteuler.net/problem=735
Möbius function: Wikipedia — Möbius function
Squarefree integer: Wikipedia — Squarefree integer
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method

Problem 735 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128_t;

constexpr u64 kDefaultN = 1'000'000'000'000ULL;

struct Options {
    u64 n = kDefaultN;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }
    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            options.n = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

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

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

// Computes sum_{m=L}^X floor(X/m).
u64 floor_sum_from(const u64 x, const u64 l) {
    if (l > x) {
        return 0ULL;
    }

    u128 total = 0U;
    if (static_cast<u128>(l) * static_cast<u128>(l) > static_cast<u128>(x)) {
        const u64 qmax = x / l;
        for (u64 q = 1ULL; q <= qmax; ++q) {
            u64 hi = x / q;
            u64 lo = x / (q + 1ULL) + 1ULL;
            if (hi < l) {
                continue;
            }
            if (lo < l) {
                lo = l;
            }
            if (lo <= hi) {
                total += static_cast<u128>(q) * static_cast<u128>(hi - lo + 1ULL);
            }
        }
    } else {
        for (u64 m = l; m <= x;) {
            const u64 q = x / m;
            const u64 r = x / q;
            total += static_cast<u128>(q) * static_cast<u128>(r - m + 1ULL);
            m = r + 1ULL;
        }
    }

    return static_cast<u64>(total);
}

struct MuDivCache {
    std::vector<std::vector<std::pair<u32, std::int8_t>>> divs;
};

MuDivCache build_mu_squarefree_divisors(const int nmax) {
    MuDivCache out;
    out.divs.resize(static_cast<std::size_t>(nmax) + 1ULL);
    out.divs[1].push_back({1U, 1});

    std::vector<int> spf(static_cast<std::size_t>(nmax) + 1ULL, 0);
    for (int i = 2; i <= nmax; ++i) {
        if (spf[static_cast<std::size_t>(i)] != 0) {
            continue;
        }
        spf[static_cast<std::size_t>(i)] = i;
        if (static_cast<u64>(i) * static_cast<u64>(i) > static_cast<u64>(nmax)) {
            continue;
        }
        for (u64 j = static_cast<u64>(i) * static_cast<u64>(i); j <= static_cast<u64>(nmax);
             j += static_cast<u64>(i)) {
            if (spf[static_cast<std::size_t>(j)] == 0) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    for (int i = 2; i <= nmax; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
        }
    }

    std::vector<int> primes;
    for (int x = 2; x <= nmax; ++x) {
        int y = x;
        primes.clear();
        while (y > 1) {
            const int p = spf[static_cast<std::size_t>(y)];
            primes.push_back(p);
            while (y % p == 0) {
                y /= p;
            }
        }

        const int k = static_cast<int>(primes.size());
        const int total_masks = 1 << k;
        auto& bucket = out.divs[static_cast<std::size_t>(x)];
        bucket.reserve(static_cast<std::size_t>(total_masks));

        for (int mask = 0; mask < total_masks; ++mask) {
            u32 d = 1U;
            int bit_count = 0;
            for (int i = 0; i < k; ++i) {
                if (((mask >> i) & 1) == 0) {
                    continue;
                }
                d *= static_cast<u32>(primes[static_cast<std::size_t>(i)]);
                ++bit_count;
            }
            const std::int8_t mu = static_cast<std::int8_t>((bit_count & 1) ? -1 : 1);
            bucket.push_back({d, mu});
        }
    }

    return out;
}

struct WorkerState {
    const MuDivCache* mu_cache = nullptr;
    u64 n = 0ULL;
    u64 odd_u_max = 0ULL;
    u64 even_w_max = 0ULL;
    std::atomic<u64> next_odd_u{1ULL};
    std::atomic<u64> next_even_w{1ULL};
};

i128 process_odd_u_block(const WorkerState& state, const u64 begin_u, const u64 end_u_exclusive) {
    i128 local = 0;
    for (u64 u = begin_u; u < end_u_exclusive; u += 2ULL) {
        if (u > state.odd_u_max) {
            break;
        }
        const auto& divisors = state.mu_cache->divs[static_cast<std::size_t>(u)];
        for (const auto& [d32, mu8] : divisors) {
            const u64 d = static_cast<u64>(d32);
            const u64 x = state.n / (u * d);
            if (x == 0ULL) {
                continue;
            }
            const u64 l = (u + d - 1ULL) / d; // ceil(u/d)
            const u64 contrib = floor_sum_from(x, l);
            if (mu8 > 0) {
                local += static_cast<i128>(contrib);
            } else {
                local -= static_cast<i128>(contrib);
            }
        }
    }
    return local;
}

i128 process_even_w_block(const WorkerState& state, const u64 begin_w, const u64 end_w_exclusive) {
    i128 local = 0;
    for (u64 w = begin_w; w < end_w_exclusive; ++w) {
        if (w > state.even_w_max) {
            break;
        }
        const auto& divisors = state.mu_cache->divs[static_cast<std::size_t>(w)];
        for (const auto& [d32, mu8] : divisors) {
            const u64 d = static_cast<u64>(d32);
            if ((d & 1ULL) == 0ULL) {
                continue;
            }

            const u64 x = state.n / (w * d);
            if (x == 0ULL) {
                continue;
            }
            const u64 l = (2ULL * w + d - 1ULL) / d; // ceil(2w/d)
            if (l > x) {
                continue;
            }

            // odd m only: sum_{odd m >= l} floor(x/m)
            const u64 odd_part =
                floor_sum_from(x, l) - floor_sum_from(x >> 1ULL, (l + 1ULL) >> 1ULL);
            if (mu8 > 0) {
                local += static_cast<i128>(odd_part);
            } else {
                local -= static_cast<i128>(odd_part);
            }
        }
    }
    return local;
}

u128 compute_f_sum(const u64 n, const MuDivCache& mu_cache, const unsigned threads) {
    const u64 odd_u_max = isqrt_u64(n);
    const u64 even_w_max = isqrt_u64(n / 2ULL);

    if (threads <= 1U) {
        const WorkerState state{&mu_cache, n, odd_u_max, even_w_max};
        i128 total = 0;
        total += process_odd_u_block(state, 1ULL, odd_u_max + 2ULL);
        total += process_even_w_block(state, 1ULL, even_w_max + 1ULL);
        return static_cast<u128>(total);
    }

    WorkerState state;
    state.mu_cache = &mu_cache;
    state.n = n;
    state.odd_u_max = odd_u_max;
    state.even_w_max = even_w_max;
    state.next_odd_u.store(1ULL, std::memory_order_relaxed);
    state.next_even_w.store(1ULL, std::memory_order_relaxed);

    constexpr u64 kOddBlockSize = 256ULL; // counts odd u's, stride 2
    constexpr u64 kEvenBlockSize = 512ULL;

    std::vector<std::thread> pool;
    std::vector<i128> partial(static_cast<std::size_t>(threads), 0);
    pool.reserve(static_cast<std::size_t>(threads));

    for (unsigned tid = 0U; tid < threads; ++tid) {
        pool.emplace_back([&, tid]() {
            i128 local = 0;
            while (true) {
                bool did_work = false;

                const u64 odd_start =
                    state.next_odd_u.fetch_add(2ULL * kOddBlockSize, std::memory_order_relaxed);
                if (odd_start <= state.odd_u_max) {
                    const u64 odd_end = odd_start + 2ULL * kOddBlockSize;
                    local += process_odd_u_block(state, odd_start, odd_end);
                    did_work = true;
                }

                const u64 even_start =
                    state.next_even_w.fetch_add(kEvenBlockSize, std::memory_order_relaxed);
                if (even_start <= state.even_w_max) {
                    const u64 even_end = even_start + kEvenBlockSize;
                    local += process_even_w_block(state, even_start, even_end);
                    did_work = true;
                }

                if (!did_work) {
                    break;
                }
            }
            partial[static_cast<std::size_t>(tid)] = local;
        });
    }

    for (std::thread& th : pool) {
        th.join();
    }

    i128 total = 0;
    for (const i128 value : partial) {
        total += value;
    }

    return static_cast<u128>(total);
}

u64 brute_f(const u64 n) {
    u64 total = 0ULL;
    for (u64 x = 1ULL; x <= n; ++x) {
        const u128 m = static_cast<u128>(2ULL) * static_cast<u128>(x) * static_cast<u128>(x);
        u64 c = 0ULL;
        for (u64 d = 1ULL; d <= x; ++d) {
            if (m % static_cast<u128>(d) == 0U) {
                ++c;
            }
        }
        total += c;
    }
    return total;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const u64 odd_u_max,
                             const u64 even_w_max) {
    if (!allow_multithreading) {
        return 1U;
    }

    const u64 workload = odd_u_max / 2ULL + even_w_max;
    if (workload < 200'000ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
        threads = std::min<unsigned>(threads, 16U);
    }
    return std::max(1U, threads);
}

bool run_checkpoints(const MuDivCache& mu_cache, const unsigned threads) {
    struct Checkpoint {
        u64 n;
        u64 expected;
    };
    const std::vector<Checkpoint> checkpoints = {
        {15ULL, 63ULL},
        {1'000ULL, 15'066ULL},
    };

    for (const Checkpoint cp : checkpoints) {
        const u128 got = compute_f_sum(cp.n, mu_cache, 1U);
        if (got != static_cast<u128>(cp.expected)) {
            std::cerr << "Checkpoint failed: F(" << cp.n << ") expected " << cp.expected
                      << ", got " << to_string_u128(got) << '\n';
            return false;
        }
    }

    const u64 brute_limit = 120ULL;
    const u64 brute_expected = brute_f(brute_limit);
    const u128 fast = compute_f_sum(brute_limit, mu_cache, 1U);
    if (fast != static_cast<u128>(brute_expected)) {
        std::cerr << "Brute-force cross-check failed at N=" << brute_limit
                  << ": expected " << brute_expected << ", got " << to_string_u128(fast) << '\n';
        return false;
    }

    if (threads > 1U) {
        const u64 thread_check_n = 10'000'000ULL;
        const u128 single = compute_f_sum(thread_check_n, mu_cache, 1U);
        const u128 multi = compute_f_sum(thread_check_n, mu_cache, threads);
        if (single != multi) {
            std::cerr << "Thread consistency failed at N=" << thread_check_n
                      << ": single=" << to_string_u128(single)
                      << ", multi=" << to_string_u128(multi) << '\n';
            return false;
        }
    }

    std::cout << "Validation checkpoints passed.\n";
    return true;
}

} // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const u64 odd_u_max = isqrt_u64(options.n);
    const u64 even_w_max = isqrt_u64(options.n / 2ULL);
    const u64 limit_u64 = std::max(odd_u_max, even_w_max) + 8ULL;
    if (limit_u64 > static_cast<u64>(std::numeric_limits<int>::max())) {
        std::cerr << "N is too large for this implementation.\n";
        return 1;
    }
    const int divisor_cache_limit = static_cast<int>(limit_u64);
    const MuDivCache mu_cache = build_mu_squarefree_divisors(divisor_cache_limit);

    const unsigned threads = choose_thread_count(
        options.allow_multithreading, options.requested_threads, odd_u_max, even_w_max);

    if (options.run_checkpoints && !run_checkpoints(mu_cache, threads)) {
        return 1;
    }

    const u128 answer = compute_f_sum(options.n, mu_cache, threads);
    std::cout << to_string_u128(answer) << '\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 Euler735 {
    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("Euler735.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(".euler735_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 Euler735 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("Euler735 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("Euler735 C++ bridge produced empty output.");
        }
        return parsed;
    }

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