Problem 410: Circle and Tangent Line

Project Euler Problem 410 Solution

EulerSolve provides an optimized solution for Project Euler Problem 410, Circle and Tangent Line, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The geometric statement can be converted into an arithmetic counting problem. For a fixed bound pair $(A,B)$, every admissible tangent configuration is encoded by a positive integer $n \le A$, a positive integer $b \le B$, and a factor pair of the square $b^2$. The direct formulation is still too slow, but it exposes enough structure to turn the geometry into a divisor sum. The official task requires two such counts, for $(10^8,10^9)$ and $(10^9,10^8)$, and adds them together. The whole challenge is therefore to derive a fast closed form for a single subproblem $F(A,B)$. Mathematical Approach Step 1: Arithmetic form of the tangent condition The direct counting model shows that $F(A,B)$ can be written in terms of positive integers satisfying $$1 \le n \le A,\qquad 1 \le b \le B,\qquad cd=b^2,$$ together with the two arithmetic constraints $$b \mid n(c+d),\qquad \frac{n(c+d)}{b}\equiv d-c \pmod 2.$$ Whenever these conditions hold, there are two symmetric sign choices, so each admissible arithmetic encoding contributes exactly $2$ geometric configurations. This is the only part where the geometry still appears; from here onward the problem is pure number theory. Step 2: Every factor pair of a square has a coprime-square core Let $\lambda=\gcd(c,d)$. Then $c=\lambda c_0$, $d=\lambda d_0$, and $\gcd(c_0,d_0)=1$....

Detailed mathematical approach

Problem Summary

The geometric statement can be converted into an arithmetic counting problem. For a fixed bound pair $(A,B)$, every admissible tangent configuration is encoded by a positive integer $n \le A$, a positive integer $b \le B$, and a factor pair of the square $b^2$. The direct formulation is still too slow, but it exposes enough structure to turn the geometry into a divisor sum.

The official task requires two such counts, for $(10^8,10^9)$ and $(10^9,10^8)$, and adds them together. The whole challenge is therefore to derive a fast closed form for a single subproblem $F(A,B)$.

Mathematical Approach

Step 1: Arithmetic form of the tangent condition

The direct counting model shows that $F(A,B)$ can be written in terms of positive integers satisfying

$$1 \le n \le A,\qquad 1 \le b \le B,\qquad cd=b^2,$$

together with the two arithmetic constraints

$$b \mid n(c+d),\qquad \frac{n(c+d)}{b}\equiv d-c \pmod 2.$$

Whenever these conditions hold, there are two symmetric sign choices, so each admissible arithmetic encoding contributes exactly $2$ geometric configurations. This is the only part where the geometry still appears; from here onward the problem is pure number theory.

Step 2: Every factor pair of a square has a coprime-square core

Let $\lambda=\gcd(c,d)$. Then $c=\lambda c_0$, $d=\lambda d_0$, and $\gcd(c_0,d_0)=1$. Since $cd=b^2$, we obtain

$$c_0d_0=\left(\frac{b}{\lambda}\right)^2.$$

A product of two coprime positive integers can be a square only if each factor is already a square. Hence there exist coprime integers $p,q \ge 1$ such that

$$c=\lambda p^2,\qquad d=\lambda q^2,\qquad \gcd(p,q)=1,$$

and therefore

$$b=\lambda pq.$$

This is the key normalization step: the square condition is no longer hidden inside divisors of $b^2$; it is explicit in the parameters $(\lambda,p,q)$.

Step 3: The divisibility condition collapses to $pq \mid n$

Substituting the decomposition above into the divisibility condition gives

$$\lambda pq \mid n\lambda(p^2+q^2),$$

so equivalently

$$pq \mid n(p^2+q^2).$$

Now $\gcd(p,q)=1$, and no prime dividing $pq$ can divide $p^2+q^2$. Therefore

$$\gcd(pq,p^2+q^2)=1,$$

which forces

$$pq \mid n.$$

Write

$$m=pq,\qquad n=mt.$$

Then $t$ can range through

$$1 \le t \le \left\lfloor\frac{A}{m}\right\rfloor,$$

while the scale factor $\lambda$ is constrained only by $b=\lambda m \le B$, hence

$$1 \le \lambda \le \left\lfloor\frac{B}{m}\right\rfloor.$$

Step 4: Why the weight is $2^{\omega(m)}$

Fix $m$. We must count ordered coprime pairs $(p,q)$ such that $pq=m$. If

$$m=\prod_{i=1}^{r}\ell_i^{e_i},$$

then coprimality means each whole prime power $\ell_i^{e_i}$ must go entirely to $p$ or entirely to $q$. Every distinct prime gives exactly two choices, so the number of ordered pairs is

$$2^r=2^{\omega(m)},$$

where $\omega(m)$ is the number of distinct prime divisors of $m$. This explains the multiplicative weight used in the implementation.

Step 5: Odd $m$ and even $m$ behave differently

After substitution, the parity condition becomes

$$\frac{n(c+d)}{b}=t(p^2+q^2),\qquad d-c=\lambda(q^2-p^2).$$

We need these two quantities to have the same parity. Since

$$p^2+q^2\equiv q^2-p^2 \pmod 2,$$

the condition is equivalent to

$$ (t-\lambda)(p^2+q^2)\equiv 0 \pmod 2.$$

If $m$ is odd, then both $p$ and $q$ are odd, so $p^2+q^2$ is even. The parity restriction disappears completely. For each of the $2^{\omega(m)}$ coprime factorizations, every pair $(t,\lambda)$ is valid, and the two sign choices contribute

$$2\left\lfloor\frac{A}{m}\right\rfloor\left\lfloor\frac{B}{m}\right\rfloor.$$

If $m$ is even, exactly one of $p,q$ is even, so $p^2+q^2$ is odd. Now we must have $t\equiv\lambda\pmod 2$. Let

$$t_{\max}=\left\lfloor\frac{A}{m}\right\rfloor,\qquad \lambda_{\max}=\left\lfloor\frac{B}{m}\right\rfloor.$$

The number of pairs $(t,\lambda)$ with equal parity is

$$\left\lceil\frac{t_{\max}}{2}\right\rceil\left\lceil\frac{\lambda_{\max}}{2}\right\rceil+\left\lfloor\frac{t_{\max}}{2}\right\rfloor\left\lfloor\frac{\lambda_{\max}}{2}\right\rfloor =\frac{t_{\max}\lambda_{\max}+(t_{\max}\bmod 2)(\lambda_{\max}\bmod 2)}{2}.$$

After multiplying by the two sign choices, the even case contributes

$$t_{\max}\lambda_{\max}+(t_{\max}\bmod 2)(\lambda_{\max}\bmod 2).$$

Step 6: Final summation formula

Define the parity-dependent unit contribution by

$$U_m(A,B)= \begin{cases} 2\left\lfloor\frac{A}{m}\right\rfloor\left\lfloor\frac{B}{m}\right\rfloor, & m\ \mathrm{odd},\\[6pt] \left\lfloor\frac{A}{m}\right\rfloor\left\lfloor\frac{B}{m}\right\rfloor+ \left(\left\lfloor\frac{A}{m}\right\rfloor\bmod 2\right)\left(\left\lfloor\frac{B}{m}\right\rfloor\bmod 2\right), & m\ \mathrm{even}. \end{cases}$$

Then the whole subproblem is

$$\boxed{F(A,B)=\sum_{m=1}^{\min(A,B)}2^{\omega(m)}\,U_m(A,B).}$$

The required total is

$$F(10^8,10^9)+F(10^9,10^8).$$

The summand is symmetric in $A$ and $B$, so the two terms are equal, but writing the result as a sum of two evaluations mirrors the original statement and the implementations.

Worked Example: $F(2,10)=52$

Only $m=1$ and $m=2$ can contribute.

For $m=1$, we have $\omega(1)=0$, so the weight is $1$. Since $1$ is odd,

$$U_1(2,10)=2\cdot 2\cdot 10=40.$$

For $m=2$, we have $\omega(2)=1$, so the weight is $2$. Since $2$ is even,

$$U_2(2,10)=\left\lfloor\frac{2}{2}\right\rfloor\left\lfloor\frac{10}{2}\right\rfloor+ \left(\left\lfloor\frac{2}{2}\right\rfloor\bmod 2\right)\left(\left\lfloor\frac{10}{2}\right\rfloor\bmod 2\right) =1\cdot 5+1\cdot 1=6.$$

Therefore

$$F(2,10)=1\cdot 40+2\cdot 6=52,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations first precompute $\omega(m)$ for every $m$ up to the largest required limit by a sieve: each prime visits all of its multiples once, so every integer accumulates its number of distinct prime divisors. Then the implementation evaluates the summation above term by term, using only floor divisions, parity checks, a bit shift for $2^{\omega(m)}$, and integer accumulation.

The same scalar routine is applied to the two official bound pairs. Because each $m$-term is independent, the outer sum is naturally parallelizable; the C++ and Java implementations exploit that independence, while the Python implementation keeps the same mathematics in a straightforward single-threaded loop.

Complexity Analysis

Let $L=\max(\min(10^8,10^9),\min(10^9,10^8))=10^8$. Building the distinct-prime-factor table costs $O(L\log\log L)$ time and $O(L)$ memory. The final summation is $O(L)$ time because every $m$ contributes a constant-time term. Thus the overall method runs in $O(L\log\log L)$ time and uses $O(L)$ space.

References

Problem page: https://projecteuler.net/problem=410
Prime omega function: Wikipedia — Prime omega function
Fundamental theorem of arithmetic: Wikipedia — Fundamental theorem of arithmetic
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes
Hardy, G. H., Wright, E. M. An Introduction to the Theory of Numbers. Sections on unique factorization and arithmetic functions.

Problem 410 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u8 = std::uint8_t;
using u128 = unsigned __int128;

constexpr u64 kDefaultR1 = 100'000'000ULL;
constexpr u64 kDefaultX1 = 1'000'000'000ULL;
constexpr u64 kDefaultR2 = 1'000'000'000ULL;
constexpr u64 kDefaultX2 = 100'000'000ULL;

constexpr u64 kCheckpointR1 = 1ULL;
constexpr u64 kCheckpointX1 = 5ULL;
constexpr u64 kCheckpointExpected1 = 10ULL;
constexpr u64 kCheckpointR2 = 2ULL;
constexpr u64 kCheckpointX2 = 10ULL;
constexpr u64 kCheckpointExpected2 = 52ULL;
constexpr u64 kCheckpointR3 = 10ULL;
constexpr u64 kCheckpointX3 = 100ULL;
constexpr u64 kCheckpointExpected3 = 3384ULL;

constexpr u64 kBruteCheckR = 12ULL;
constexpr u64 kBruteCheckX = 20ULL;
constexpr u64 kBruteCheckExpected = 824ULL;

constexpr u64 kThreadConsistencyR = 200'000ULL;
constexpr u64 kThreadConsistencyX = 300'000ULL;

struct Options {
    u64 r1 = kDefaultR1;
    u64 x1 = kDefaultX1;
    u64 r2 = kDefaultR2;
    u64 x2 = kDefaultX2;
    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) != 0) {
        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(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, "--r1=", parsed_u64)) {
            options.r1 = parsed_u64;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--x1=", parsed_u64)) {
            options.x1 = parsed_u64;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--r2=", parsed_u64)) {
            options.r2 = parsed_u64;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--x2=", parsed_u64)) {
            options.x2 = 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;
    }

    if (options.r1 == 0ULL || options.x1 == 0ULL || options.r2 == 0ULL || options.x2 == 0ULL) {
        std::cerr << "All bounds must be at least 1.\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }

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

std::vector<u8> build_distinct_prime_counts(u64 limit) {
    std::vector<u8> omega(static_cast<std::size_t>(limit) + 1ULL, 0U);
    if (limit < 2ULL) {
        return omega;
    }

    std::vector<bool> composite((static_cast<std::size_t>(limit) >> 1ULL) + 1ULL, false);
    const u64 root = static_cast<u64>(std::sqrt(static_cast<long double>(limit)));

    for (u64 p = 3ULL; p <= root; p += 2ULL) {
        if (composite[static_cast<std::size_t>(p >> 1ULL)]) {
            continue;
        }

        const u64 step = p << 1ULL;
        for (u64 j = p * p; j <= limit; j += step) {
            composite[static_cast<std::size_t>(j >> 1ULL)] = true;
        }
    }

    for (u64 m = 2ULL; m <= limit; m += 2ULL) {
        ++omega[static_cast<std::size_t>(m)];
    }
    for (u64 p = 3ULL; p <= limit; p += 2ULL) {
        if (composite[static_cast<std::size_t>(p >> 1ULL)]) {
            continue;
        }
        for (u64 m = p; m <= limit; m += p) {
            ++omega[static_cast<std::size_t>(m)];
        }
    }

    return omega;
}

u64 brute_force_reference(u64 R, u64 X) {
    u64 total = 0ULL;

    for (u64 r = 1ULL; r <= R; ++r) {
        for (u64 a = 1ULL; a <= X; ++a) {
            const u64 a2 = a * a;

            for (u64 u = 1ULL; u <= a2; ++u) {
                if (a2 % u != 0ULL) {
                    continue;
                }
                const u64 v = a2 / u;
                const long long d = static_cast<long long>(v) - static_cast<long long>(u);
                const u64 k = u + v;

                const u64 rk = r * k;
                if (rk % a != 0ULL) {
                    continue;
                }

                const long long s_abs = static_cast<long long>(rk / a);
                if (((s_abs - d) & 1LL) != 0LL) {
                    continue;
                }

                total += 2ULL;  // s = +s_abs and s = -s_abs
            }
        }
    }

    return total;
}

u128 solve_F_with_omega(u64 R,
                        u64 X,
                        const std::vector<u8>& omega,
                        bool allow_multithreading,
                        unsigned requested_threads) {
    const u64 limit = std::min(R, X);
    if (limit == 0ULL) {
        return 0;
    }

    const unsigned threads = choose_thread_count(
        allow_multithreading, requested_threads, static_cast<std::size_t>(limit));
    std::vector<u128> partial(threads, 0);

    auto worker = [&](unsigned tid) {
        u128 local = 0;

        for (u64 m = 1ULL + static_cast<u64>(tid); m <= limit; m += threads) {
            const u64 t = R / m;
            const u64 g = X / m;
            const u64 weight = 1ULL << omega[static_cast<std::size_t>(m)];

            u64 unit = 0ULL;
            if ((m & 1ULL) != 0ULL) {
                unit = 2ULL * t * g;
            } else {
                unit = t * g + ((t & 1ULL) & (g & 1ULL));
            }

            local += static_cast<u128>(weight) * static_cast<u128>(unit);
        }

        partial[tid] = local;
    };

    std::vector<std::thread> pool;
    pool.reserve(threads > 0U ? threads - 1U : 0U);
    for (unsigned t = 1U; t < threads; ++t) {
        pool.emplace_back(worker, t);
    }
    worker(0U);
    for (std::thread& thread : pool) {
        thread.join();
    }

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

u128 solve_F(u64 R,
             u64 X,
             bool allow_multithreading,
             unsigned requested_threads) {
    const u64 limit = std::min(R, X);
    const std::vector<u8> omega = build_distinct_prime_counts(limit);
    return solve_F_with_omega(R, X, omega, allow_multithreading, requested_threads);
}

bool run_checkpoints(const Options& options) {
    {
        const u128 got = solve_F(kCheckpointR1,
                                 kCheckpointX1,
                                 options.allow_multithreading,
                                 options.requested_threads);
        if (got != static_cast<u128>(kCheckpointExpected1)) {
            std::cerr << "Checkpoint failed: F(" << kCheckpointR1 << ", " << kCheckpointX1
                      << ") expected " << kCheckpointExpected1
                      << ", got " << to_string_u128(got) << '\n';
            return false;
        }
        std::cout << "Checkpoint passed: F(" << kCheckpointR1 << ", " << kCheckpointX1
                  << ") = " << kCheckpointExpected1 << '\n';
    }

    {
        const u128 got = solve_F(kCheckpointR2,
                                 kCheckpointX2,
                                 options.allow_multithreading,
                                 options.requested_threads);
        if (got != static_cast<u128>(kCheckpointExpected2)) {
            std::cerr << "Checkpoint failed: F(" << kCheckpointR2 << ", " << kCheckpointX2
                      << ") expected " << kCheckpointExpected2
                      << ", got " << to_string_u128(got) << '\n';
            return false;
        }
        std::cout << "Checkpoint passed: F(" << kCheckpointR2 << ", " << kCheckpointX2
                  << ") = " << kCheckpointExpected2 << '\n';
    }

    {
        const u128 got = solve_F(kCheckpointR3,
                                 kCheckpointX3,
                                 options.allow_multithreading,
                                 options.requested_threads);
        if (got != static_cast<u128>(kCheckpointExpected3)) {
            std::cerr << "Checkpoint failed: F(" << kCheckpointR3 << ", " << kCheckpointX3
                      << ") expected " << kCheckpointExpected3
                      << ", got " << to_string_u128(got) << '\n';
            return false;
        }
        std::cout << "Checkpoint passed: F(" << kCheckpointR3 << ", " << kCheckpointX3
                  << ") = " << kCheckpointExpected3 << '\n';
    }

    {
        const u64 brute = brute_force_reference(kBruteCheckR, kBruteCheckX);
        if (brute != kBruteCheckExpected) {
            std::cerr << "Internal brute checkpoint mismatch: expected "
                      << kBruteCheckExpected << ", got " << brute << '\n';
            return false;
        }

        const u128 fast = solve_F(kBruteCheckR,
                                  kBruteCheckX,
                                  options.allow_multithreading,
                                  options.requested_threads);
        if (fast != static_cast<u128>(brute)) {
            std::cerr << "Fast-vs-brute checkpoint failed at F(" << kBruteCheckR << ", "
                      << kBruteCheckX << "): brute=" << brute
                      << ", fast=" << to_string_u128(fast) << '\n';
            return false;
        }
        std::cout << "Checkpoint passed: F(" << kBruteCheckR << ", " << kBruteCheckX
                  << ") = " << brute << " (matches brute force)\n";
    }

    if (options.allow_multithreading) {
        const u64 limit = std::min(kThreadConsistencyR, kThreadConsistencyX);
        const std::vector<u8> omega = build_distinct_prime_counts(limit);

        const u128 multi = solve_F_with_omega(kThreadConsistencyR,
                                              kThreadConsistencyX,
                                              omega,
                                              true,
                                              options.requested_threads);
        const u128 single = solve_F_with_omega(kThreadConsistencyR,
                                               kThreadConsistencyX,
                                               omega,
                                               false,
                                               options.requested_threads);
        if (multi != single) {
            std::cerr << "Thread consistency failed at F(" << kThreadConsistencyR << ", "
                      << kThreadConsistencyX << "): single=" << to_string_u128(single)
                      << ", multi=" << to_string_u128(multi) << '\n';
            return false;
        }
        std::cout << "Checkpoint passed: thread consistency at F(" << kThreadConsistencyR
                  << ", " << kThreadConsistencyX << ")\n";
    }

    return true;
}

}  // namespace

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

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    try {
        const auto start_time = std::chrono::steady_clock::now();

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

        const u64 limit1 = std::min(options.r1, options.x1);
        const u64 limit2 = std::min(options.r2, options.x2);
        const u64 max_limit = std::max(limit1, limit2);
        const std::vector<u8> omega = build_distinct_prime_counts(max_limit);

        const u128 f1 = solve_F_with_omega(options.r1,
                                           options.x1,
                                           omega,
                                           options.allow_multithreading,
                                           options.requested_threads);

        u128 f2 = 0;
        if (options.r2 == options.x1 && options.x2 == options.r1) {
            f2 = f1;
        } else {
            f2 = solve_F_with_omega(options.r2,
                                    options.x2,
                                    omega,
                                    options.allow_multithreading,
                                    options.requested_threads);
        }

        const u128 answer = f1 + f2;

        const auto end_time = std::chrono::steady_clock::now();
        const std::chrono::duration<double> elapsed = end_time - start_time;

        std::cout << "F(" << options.r1 << ", " << options.x1 << ") = "
                  << to_string_u128(f1) << '\n';
        std::cout << "F(" << options.r2 << ", " << options.x2 << ") = "
                  << to_string_u128(f2) << '\n';
        std::cout << "Answer = " << to_string_u128(answer) << '\n';
        std::cout << "Elapsed: " << elapsed.count() << " seconds\n";
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

def solve():
    R1, X1 = 100000000, 1000000000
    R2, X2 = 1000000000, 100000000
    limit = max(min(R1, X1), min(R2, X2))

    # Build distinct prime factor count (omega) via sieve
    omega = bytearray(limit + 1)
    is_comp = bytearray((limit >> 1) + 1)
    root = int(limit**0.5)
    for p in range(3, root + 1, 2):
        if is_comp[p >> 1]: continue
        for j in range(p * p, limit + 1, p << 1):
            is_comp[j >> 1] = 1
    for m in range(2, limit + 1, 2):
        omega[m] += 1
    for p in range(3, limit + 1, 2):
        if is_comp[p >> 1]: continue
        for m in range(p, limit + 1, p):
            omega[m] += 1

    def solve_F(R, X):
        lim = min(R, X)
        total = 0
        for m in range(1, lim + 1):
            t = R // m
            g = X // m
            w = 1 << omega[m]
            if m & 1:
                unit = 2 * t * g
            else:
                unit = t * g + ((t & 1) & (g & 1))
            total += w * unit
        return total

    return str(solve_F(R1, X1) + solve_F(R2, X2))

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

Java

import java.math.BigInteger;
import java.util.stream.IntStream;

public class Euler410 {
    static byte[] buildDistinctPrimeCounts(int limit) {
        byte[] omega = new byte[limit + 1];
        if (limit < 2)
            return omega;

        byte[] composite = new byte[(limit >> 1) + 1];
        int root = (int) Math.sqrt(limit);

        for (int p = 3; p <= root; p += 2) {
            if (composite[p >> 1] != 0)
                continue;
            int step = p << 1;
            for (int j = p * p; j <= limit; j += step) {
                composite[j >> 1] = 1;
            }
        }

        for (int m = 2; m <= limit; m += 2) {
            omega[m]++;
        }
        for (int p = 3; p <= limit; p += 2) {
            if (composite[p >> 1] != 0)
                continue;
            for (int m = p; m <= limit; m += p) {
                omega[m]++;
            }
        }

        return omega;
    }

    static BigInteger solveF(long R, long X, byte[] omega) {
        int limit = (int) Math.min(R, X);
        if (limit == 0)
            return BigInteger.ZERO;

        int numChunks = Runtime.getRuntime().availableProcessors();
        long chunkSize = (limit + numChunks - 1) / numChunks;

        BigInteger total = IntStream.range(0, numChunks)
                .parallel()
                .mapToObj(i -> {
                    int start = (int) (i * chunkSize + 1);
                    int end = (int) Math.min((i + 1) * chunkSize, limit);

                    long sumHi = 0;
                    long sumLo = 0;

                    for (int m = start; m <= end; m++) {
                        long t = R / m;
                        long g = X / m;
                        long weight = 1L << omega[m];
                        long unit;
                        if ((m & 1) != 0) {
                            unit = 2L * t * g;
                        } else {
                            unit = t * g + ((t & 1) & (g & 1));
                        }

                        long term = weight * unit;
                        sumLo += term;
                        if (Long.compareUnsigned(sumLo, term) < 0) {
                            sumHi++;
                        }
                    }

                    BigInteger hi = BigInteger.valueOf(sumHi).shiftLeft(64);
                    BigInteger lo = BigInteger.valueOf(sumLo);
                    if (sumLo < 0) {
                        lo = lo.add(BigInteger.ONE.shiftLeft(64));
                    }
                    return hi.add(lo);
                })
                .reduce(BigInteger.ZERO, BigInteger::add);

        return total;
    }

    static String solve() {
        long r1 = 100000000L;
        long x1 = 1000000000L;
        long r2 = 1000000000L;
        long x2 = 100000000L;

        int maxLimit = (int) Math.max(Math.min(r1, x1), Math.min(r2, x2));
        byte[] omega = buildDistinctPrimeCounts(maxLimit);

        BigInteger f1 = solveF(r1, x1, omega);
        BigInteger f2;
        if (r2 == x1 && x2 == r1) {
            f2 = f1;
        } else {
            f2 = solveF(r2, x2, omega);
        }

        return f1.add(f2).toString();
    }

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