Problem 233: Lattice Points on a Circle

Project Euler Problem 233 Solution

EulerSolve provides an optimized solution for Project Euler Problem 233, Lattice Points on a Circle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $N$, the problem considers the circle through $(0,0)$, $(N,0)$, $(0,N)$, and $(N,N)$, and defines $f(N)$ as the number of lattice points on that circle. The goal is to compute $$S(10^{11})=\sum_{\substack{N\le 10^{11}\\f(N)=420}} N.$$ A direct scan up to $10^{11}$ is impossible. The successful route is to convert the geometry into a sum-of-two-squares counting problem, extract an exact multiplicative formula for $f(N)$, and then classify all admissible prime-exponent patterns. Mathematical Approach The central observation is that the original circle can be rewritten so that its lattice points are counted by representations of $N^2$ as a sum of two squares. From the original circle to $a^2+b^2=N^2$ The circle through the four corners of the $N\times N$ square has center $(N/2,N/2)$ and radius $N/\sqrt2$, so its equation is $$\left(x-\frac N2\right)^2+\left(y-\frac N2\right)^2=\frac{N^2}{2}.$$ Introduce the linear change of variables $$a=x+y-N,\qquad b=x-y.$$ Substituting and simplifying gives $$a^2+b^2=N^2.$$ This transformation is bijective on lattice points: from $a$ and $b$ we recover $$x=\frac{N+a+b}{2},\qquad y=\frac{N+a-b}{2},$$ and the parity condition is automatically satisfied whenever $a^2+b^2=N^2$. Therefore $f(N)$ is exactly the number of integer pairs $(a,b)$ satisfying that equation....

Detailed mathematical approach

Problem Summary

For each positive integer $N$, the problem considers the circle through $(0,0)$, $(N,0)$, $(0,N)$, and $(N,N)$, and defines $f(N)$ as the number of lattice points on that circle. The goal is to compute

$$S(10^{11})=\sum_{\substack{N\le 10^{11}\\f(N)=420}} N.$$

A direct scan up to $10^{11}$ is impossible. The successful route is to convert the geometry into a sum-of-two-squares counting problem, extract an exact multiplicative formula for $f(N)$, and then classify all admissible prime-exponent patterns.

Mathematical Approach

The central observation is that the original circle can be rewritten so that its lattice points are counted by representations of $N^2$ as a sum of two squares.

From the original circle to $a^2+b^2=N^2$

The circle through the four corners of the $N\times N$ square has center $(N/2,N/2)$ and radius $N/\sqrt2$, so its equation is

$$\left(x-\frac N2\right)^2+\left(y-\frac N2\right)^2=\frac{N^2}{2}.$$

Introduce the linear change of variables

$$a=x+y-N,\qquad b=x-y.$$

Substituting and simplifying gives

$$a^2+b^2=N^2.$$

This transformation is bijective on lattice points: from $a$ and $b$ we recover

$$x=\frac{N+a+b}{2},\qquad y=\frac{N+a-b}{2},$$

and the parity condition is automatically satisfied whenever $a^2+b^2=N^2$. Therefore $f(N)$ is exactly the number of integer pairs $(a,b)$ satisfying that equation. In standard notation,

$$f(N)=r_2(N^2),$$

where $r_2(m)$ counts representations of $m$ as a sum of two squares with order and sign included.

Prime factorization and the formula for $f(N)$

Write the factorization of $N$ as

$$N=2^\alpha\prod_i p_i^{e_i}\prod_j q_j^{f_j},$$

where every $p_i\equiv 1\pmod 4$ and every $q_j\equiv 3\pmod 4$. Then

$$N^2=2^{2\alpha}\prod_i p_i^{2e_i}\prod_j q_j^{2f_j}.$$

The sum-of-two-squares theorem implies that primes $q_j\equiv 3\pmod 4$ do not affect the count here, because they appear with even exponents in $N^2$. Only the $1\bmod 4$ primes contribute, and the count becomes

$$f(N)=r_2(N^2)=4\prod_i(2e_i+1).$$

This is the exact arithmetic object computed by the implementations, and it is also why the geometry can be cross-checked by factoring $N$ instead of drawing the circle.

Solving $f(N)=420$

The target condition is

$$4\prod_i(2e_i+1)=420,$$

so we must solve

$$\prod_i(2e_i+1)=105=3\cdot5\cdot7.$$

Each factor $2e_i+1$ is an odd integer at least $3$. That leaves only five possible exponent patterns:

$$105 \Rightarrow e=52,$$

$$35\cdot3 \Rightarrow (e_1,e_2)=(17,1),$$

$$21\cdot5 \Rightarrow (e_1,e_2)=(10,2),$$

$$15\cdot7 \Rightarrow (e_1,e_2)=(7,3),$$

$$3\cdot5\cdot7 \Rightarrow (e_1,e_2,e_3)=(1,2,3).$$

So any valid $N$ must have its $1\bmod4$ prime exponents equal to one permutation of

$$[52],\ [17,1],\ [10,2],\ [7,3],\ [3,2,1].$$

Core numbers and allowed multipliers

Separate every valid $N$ into

$$N=c\cdot m.$$

The core $c$ contains all primes $p\equiv1\pmod4$ with one of the admissible exponent patterns above. The multiplier $m$ contains only the prime $2$ and primes $q\equiv3\pmod4$. This decomposition is unique because the two prime families are disjoint.

The crucial invariance is that

$$f(c\,m)=f(c)$$

whenever every prime factor of $m$ is $2$ or $3\bmod4$. So once a core is fixed, all allowed multipliers below the remaining limit generate further valid values of $N$.

If $\mathcal M(X)$ denotes the set of allowed multipliers not exceeding $X$, then

$$A(X)=\sum_{m\in\mathcal M(X)} m.$$

then the required sum is

$$S(10^{11})=\sum_{c\in\mathcal C} c\cdot A\!\left(\left\lfloor\frac{10^{11}}{c}\right\rfloor\right),$$

where $\mathcal C$ is the set of distinct core numbers generated from the five exponent patterns.

Worked example

Consider

$$N=5^3\cdot13^2\cdot17=359125.$$

All three primes are $1\bmod4$, and their exponents are $3$, $2$, and $1$. Therefore

$$f(N)=4(2\cdot3+1)(2\cdot2+1)(2\cdot1+1)=4\cdot7\cdot5\cdot3=420.$$

Now multiply by $2^4\cdot3\cdot7$. The new number

$$N'=2^4\cdot3\cdot7\cdot5^3\cdot13^2\cdot17$$

has the same exponents on its $1\bmod4$ primes, so it still satisfies $f(N')=420$. This is exactly the “core plus neutral multiplier” structure exploited by the algorithm.

How the Code Works

Enumerating admissible cores

The C++, Python, and Java implementations first sieve the primes $p\equiv1\pmod4$ up to a safe bound, then generate every core number compatible with the five exponent patterns. Exponents are assigned to strictly increasing $1\bmod4$ primes so that each factorization is produced once. If a partial product already exceeds the limit, or even the smallest possible completion would exceed the limit, that branch is pruned immediately.

Building the allowed-multiplier prefix sums

Once the smallest core is known, the maximum useful multiplier is $\lfloor 10^{11}/c_{\min}\rfloor$. The implementation builds a smallest-prime-factor table up to that bound and marks exactly those integers whose prime divisors are all $2$ or $3\bmod4$. A prefix array then stores $A(X)$, the sum of all allowed multipliers up to $X$, for every relevant $X$.

Adding the contributions

For each core $c$, the algorithm computes $X=\lfloor 10^{11}/c\rfloor$ and adds $c\cdot A(X)$. That single lookup accounts for every valid number of the form $c\,m$ without iterating over all such multiples individually. The C++ and Java implementations split this final sweep across several worker threads, while the Python implementation performs the same summation either serially or with process-level parallel work.

The C++ implementation also includes three consistency checks: it verifies the sample value $f(10000)=36$, compares the arithmetic formula against direct geometric counting for small $N$, and confirms that the fast summation matches a brute-force summation on a smaller limit.

Complexity Analysis

Let $C$ be the number of distinct core numbers and let $M=\lfloor 10^{11}/c_{\min}\rfloor$, where $c_{\min}=5^3\cdot13^2\cdot17$ is the smallest valid core. Core generation is very small compared with the original search space, because only five exponent patterns exist and the recursive search is heavily pruned by monotonic growth of prime powers.

The sieve, validity table, and prefix sums over multipliers cost $O(M)$ time and $O(M)$ memory. The final accumulation over cores costs $O(C)$. The method therefore replaces an impossible $O(10^{11})$ search with a modest prime sieve, a compact enumeration of core numbers, and a linear preprocessing pass over the multiplier range.

Footnotes and References

Problem page: Project Euler 233 - Lattice Points on a Circle
Sum of two squares theorem: Fermat's theorem on sums of two squares
Representation count $r_2(n)$: Sum of two squares function
Gaussian integers: Gaussian integer
Prime sieving background: Sieve of Eratosthenes

Problem 233 source code

C++

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

namespace {

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

constexpr u64 kDefaultLimit = 100'000'000'000ULL;
constexpr u64 kSampleN = 10'000ULL;
constexpr u64 kSampleExpectedF = 36ULL;
constexpr u64 kGeometryCrossCheckMaxN = 400ULL;
constexpr u64 kBruteSumCrossCheckLimit = 2'000'000ULL;
constexpr u64 kThreadConsistencyLimit = 10'000'000'000ULL;
constexpr u64 kExponentOneMinOtherProduct = 21'125ULL;  // 5^3 * 13^2

struct Options {
    u64 limit = kDefaultLimit;
    bool allow_multithreading = true;
    bool run_checkpoints = true;
    unsigned requested_threads = 0;
};

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 = 0;
    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 = 0;
    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 == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }

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

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

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

    if (options.limit == 0ULL) {
        std::cerr << "--limit must be >= 1.\n";
        return false;
    }

    return true;
}

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

u64 pow_with_limit(u64 base, int exponent, u64 limit) {
    u64 result = 1ULL;
    for (int i = 0; i < exponent; ++i) {
        if (result > limit / base) {
            return limit + 1ULL;
        }
        result *= base;
    }
    return result;
}

std::vector<int> sieve_primes(int limit) {
    if (limit < 2) {
        return {};
    }

    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit) + 1ULL, 1U);
    is_prime[0] = 0U;
    is_prime[1] = 0U;

    const int root = static_cast<int>(std::sqrt(static_cast<long double>(limit)));
    for (int p = 2; p <= root; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        for (std::int64_t m = static_cast<std::int64_t>(p) * p; m <= limit; m += p) {
            is_prime[static_cast<std::size_t>(m)] = 0U;
        }
    }

    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(
        static_cast<long double>(limit) /
        std::max(1.0L, std::log(static_cast<long double>(limit)))));
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] != 0U) {
            primes.push_back(p);
        }
    }
    return primes;
}

std::vector<int> primes_1mod4_up_to(int limit) {
    const std::vector<int> primes = sieve_primes(limit);
    std::vector<int> out;
    out.reserve(primes.size() / 2ULL + 1ULL);
    for (const int p : primes) {
        if ((p % 4) == 1) {
            out.push_back(p);
        }
    }
    return out;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2'000ULL) {
        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)));
}

u64 f_value_from_factorization(u64 n, const std::vector<int>& primes) {
    u64 m = n;
    u64 product = 1ULL;

    for (const int p_int : primes) {
        const u64 p = static_cast<u64>(p_int);
        if (p > m / p) {
            break;
        }
        if ((m % p) != 0ULL) {
            continue;
        }

        int exponent = 0;
        while ((m % p) == 0ULL) {
            m /= p;
            ++exponent;
        }
        if ((p & 3ULL) == 1ULL) {
            product *= static_cast<u64>(2 * exponent + 1);
        }
    }

    if (m > 1ULL && (m & 3ULL) == 1ULL) {
        product *= 3ULL;
    }

    return 4ULL * product;
}

u64 brute_circle_count(u64 n) {
    const u64 target = n * n;
    const std::int64_t nn = static_cast<std::int64_t>(n);

    u64 count = 0ULL;
    for (std::int64_t x = -nn; x <= nn; ++x) {
        const std::int64_t rem = static_cast<std::int64_t>(target) - x * x;
        if (rem < 0) {
            continue;
        }

        const u64 y = isqrt_u64(static_cast<u64>(rem));
        if (y * y != static_cast<u64>(rem)) {
            continue;
        }
        count += (y == 0ULL ? 1ULL : 2ULL);
    }

    return count;
}

std::vector<u32> build_spf(u64 max_n) {
    std::vector<u32> spf(static_cast<std::size_t>(max_n) + 1ULL, 0U);
    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(
        static_cast<long double>(max_n) /
        std::max(1.0L, std::log(std::max(2.0L, static_cast<long double>(max_n))))));

    for (u64 i = 2ULL; i <= max_n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0U) {
            spf[static_cast<std::size_t>(i)] = static_cast<u32>(i);
            primes.push_back(static_cast<u32>(i));
        }
        for (const u32 p : primes) {
            const u64 m = static_cast<u64>(p) * i;
            if (m > max_n) {
                break;
            }
            spf[static_cast<std::size_t>(m)] = p;
            if (p == spf[static_cast<std::size_t>(i)]) {
                break;
            }
        }
    }

    return spf;
}

u128 brute_sum_direct(u64 limit) {
    if (limit == 0ULL) {
        return 0;
    }

    const std::vector<u32> spf = build_spf(limit);
    u128 sum = 0;

    for (u64 n = 1ULL; n <= limit; ++n) {
        u64 m = n;
        u64 product = 1ULL;
        while (m > 1ULL) {
            const u32 p = spf[static_cast<std::size_t>(m)];
            int exponent = 0;
            while ((m % p) == 0ULL) {
                m /= p;
                ++exponent;
            }
            if ((p & 3U) == 1U) {
                product *= static_cast<u64>(2 * exponent + 1);
            }
        }

        const u64 f = 4ULL * product;
        if (f == 420ULL) {
            sum += n;
        }
    }

    return sum;
}

bool tail_fits(const std::vector<int>& ordered_exponents,
               std::size_t next_pos,
               std::size_t next_prime_idx,
               u64 remaining_limit,
               const std::vector<int>& primes_1mod4) {
    std::size_t idx = next_prime_idx;
    for (std::size_t pos = next_pos; pos < ordered_exponents.size(); ++pos, ++idx) {
        if (idx >= primes_1mod4.size()) {
            return false;
        }
        const u64 p = static_cast<u64>(primes_1mod4[idx]);
        const u64 pe = pow_with_limit(p, ordered_exponents[pos], remaining_limit);
        if (pe > remaining_limit) {
            return false;
        }
        remaining_limit /= pe;
    }
    return true;
}

void enumerate_core_for_order(const std::vector<int>& ordered_exponents,
                              std::size_t pos,
                              std::size_t start_prime_idx,
                              u64 current_product,
                              u64 limit,
                              const std::vector<int>& primes_1mod4,
                              std::vector<u64>& out) {
    if (pos == ordered_exponents.size()) {
        out.push_back(current_product);
        return;
    }

    const std::size_t remaining_slots = ordered_exponents.size() - pos;
    if (start_prime_idx + remaining_slots > primes_1mod4.size()) {
        return;
    }

    const u64 head_limit = limit / current_product;
    for (std::size_t i = start_prime_idx; i + remaining_slots <= primes_1mod4.size(); ++i) {
        const u64 p = static_cast<u64>(primes_1mod4[i]);
        const u64 pe = pow_with_limit(p, ordered_exponents[pos], head_limit);
        if (pe > head_limit) {
            break;
        }

        const u64 next_product = current_product * pe;
        if (pos + 1ULL < ordered_exponents.size()) {
            if (!tail_fits(ordered_exponents,
                           pos + 1ULL,
                           i + 1ULL,
                           limit / next_product,
                           primes_1mod4)) {
                break;
            }
        }

        enumerate_core_for_order(ordered_exponents,
                                 pos + 1ULL,
                                 i + 1ULL,
                                 next_product,
                                 limit,
                                 primes_1mod4,
                                 out);
    }
}

std::vector<u64> enumerate_core_numbers(u64 limit, const std::vector<int>& primes_1mod4) {
    // 4 * Π(2a_i + 1) = 420  =>  Π(2a_i + 1) = 105.
    const std::vector<std::vector<int>> exponent_patterns = {
        {52},
        {17, 1},
        {10, 2},
        {7, 3},
        {3, 2, 1},
    };

    std::vector<u64> core_numbers;
    core_numbers.reserve(200'000);

    for (const auto& pattern : exponent_patterns) {
        std::vector<int> permutation = pattern;
        std::sort(permutation.begin(), permutation.end());
        do {
            enumerate_core_for_order(permutation,
                                     0,
                                     0,
                                     1ULL,
                                     limit,
                                     primes_1mod4,
                                     core_numbers);
        } while (std::next_permutation(permutation.begin(), permutation.end()));
    }

    std::sort(core_numbers.begin(), core_numbers.end());
    core_numbers.erase(std::unique(core_numbers.begin(), core_numbers.end()), core_numbers.end());
    return core_numbers;
}

std::vector<u128> build_allowed_multiplier_prefix(u64 max_multiplier) {
    std::vector<u128> prefix(static_cast<std::size_t>(max_multiplier) + 1ULL, 0);
    if (max_multiplier == 0ULL) {
        return prefix;
    }

    const std::vector<u32> spf = build_spf(max_multiplier);
    std::vector<std::uint8_t> valid(static_cast<std::size_t>(max_multiplier) + 1ULL, 0U);
    valid[1] = 1U;

    for (u64 n = 2ULL; n <= max_multiplier; ++n) {
        const u32 p = spf[static_cast<std::size_t>(n)];
        const u64 m = n / static_cast<u64>(p);
        const bool allowed_prime = (p == 2U) || ((p & 3U) == 3U);
        valid[static_cast<std::size_t>(n)] =
            static_cast<std::uint8_t>(allowed_prime && valid[static_cast<std::size_t>(m)] != 0U);
    }

    u128 running = 0;
    for (u64 n = 1ULL; n <= max_multiplier; ++n) {
        if (valid[static_cast<std::size_t>(n)] != 0U) {
            running += n;
        }
        prefix[static_cast<std::size_t>(n)] = running;
    }
    return prefix;
}

u128 accumulate_range(const std::vector<u64>& core_numbers,
                     std::size_t begin,
                     std::size_t end,
                     u64 limit,
                     const std::vector<u128>& allowed_prefix) {
    u128 local = 0;
    for (std::size_t i = begin; i < end; ++i) {
        const u64 core = core_numbers[i];
        const u64 max_multiplier = limit / core;
        local += static_cast<u128>(core) * allowed_prefix[static_cast<std::size_t>(max_multiplier)];
    }
    return local;
}

u64 estimate_prime_bound(u64 limit) {
    const u64 by_exp_one = limit / kExponentOneMinOtherProduct + 100ULL;
    return std::max<u64>(200ULL, by_exp_one);
}

u128 solve_sum(u64 limit, bool allow_multithreading, unsigned requested_threads) {
    const u64 prime_bound_u64 = estimate_prime_bound(limit);
    if (prime_bound_u64 > static_cast<u64>(std::numeric_limits<int>::max())) {
        std::cerr << "Prime bound is too large for this build.\n";
        return 0;
    }

    const std::vector<int> primes_1mod4 = primes_1mod4_up_to(static_cast<int>(prime_bound_u64));
    const std::vector<u64> core_numbers = enumerate_core_numbers(limit, primes_1mod4);
    if (core_numbers.empty()) {
        return 0;
    }

    const u64 min_core = core_numbers.front();
    const u64 max_multiplier = limit / min_core;
    const std::vector<u128> allowed_prefix = build_allowed_multiplier_prefix(max_multiplier);

    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, core_numbers.size());
    if (threads == 1U) {
        return accumulate_range(core_numbers, 0, core_numbers.size(), limit, allowed_prefix);
    }

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

    for (unsigned t = 0; t < threads; ++t) {
        const std::size_t begin = core_numbers.size() * t / threads;
        const std::size_t end = core_numbers.size() * (t + 1ULL) / threads;
        pool.emplace_back([&, t, begin, end]() {
            partial[t] = accumulate_range(core_numbers, begin, end, limit, allowed_prefix);
        });
    }
    for (std::thread& th : pool) {
        th.join();
    }

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

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

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

bool run_checkpoints(const Options& options) {
    const std::vector<int> factor_primes = sieve_primes(1'000'000);

    const u64 sample_f = f_value_from_factorization(kSampleN, factor_primes);
    if (sample_f != kSampleExpectedF) {
        std::cerr << "Checkpoint failed: f(" << kSampleN << ") expected " << kSampleExpectedF
                  << ", got " << sample_f << '\n';
        return false;
    }
    std::cout << "Checkpoint OK: f(" << kSampleN << ") = " << sample_f << '\n';

    for (u64 n = 1ULL; n <= kGeometryCrossCheckMaxN; ++n) {
        const u64 brute = brute_circle_count(n);
        const u64 formula = f_value_from_factorization(n, factor_primes);
        if (brute != formula) {
            std::cerr << "Checkpoint failed: geometry/formula mismatch at N=" << n
                      << ", brute=" << brute << ", formula=" << formula << '\n';
            return false;
        }
    }
    std::cout << "Checkpoint OK: geometry cross-check for N <= " << kGeometryCrossCheckMaxN
              << '\n';

    const u64 brute_limit = std::min<u64>(kBruteSumCrossCheckLimit, options.limit);
    const u128 brute_expected = brute_sum_direct(brute_limit);
    const u128 brute_fast = solve_sum(brute_limit, false, 1U);
    if (brute_expected != brute_fast) {
        std::cerr << "Checkpoint failed: brute sum mismatch at limit=" << brute_limit
                  << ", brute=" << to_string_u128(brute_expected)
                  << ", fast=" << to_string_u128(brute_fast) << '\n';
        return false;
    }
    std::cout << "Checkpoint OK: brute sum cross-check <= " << brute_limit
              << " gives " << to_string_u128(brute_fast) << '\n';

    if (options.allow_multithreading) {
        const u64 tc_limit = std::min<u64>(kThreadConsistencyLimit, options.limit);
        const u128 single = solve_sum(tc_limit, false, 1U);
        const u128 multi = solve_sum(tc_limit, true, options.requested_threads);
        if (single != multi) {
            std::cerr << "Checkpoint failed: thread consistency mismatch at limit=" << tc_limit
                      << ", single=" << to_string_u128(single)
                      << ", multi=" << to_string_u128(multi) << '\n';
            return false;
        }
        std::cout << "Checkpoint OK: threaded consistency at limit=" << tc_limit
                  << " gives " << to_string_u128(single) << '\n';
    }

    return true;
}

}  // namespace

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

    const auto start = std::chrono::steady_clock::now();

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

    const u128 answer =
        solve_sum(options.limit, options.allow_multithreading, options.requested_threads);

    const auto finish = std::chrono::steady_clock::now();
    const std::chrono::duration<long double> elapsed = finish - start;

    std::cout << "Answer: " << to_string_u128(answer) << '\n';
    std::cout << "S(" << options.limit << ") = " << to_string_u128(answer) << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " s\n";

    return 0;
}

Python

import math
import multiprocessing
import itertools

def isqrt(x):
    if x == 0: return 0
    return int(math.isqrt(x))

def pow_with_limit(base, exponent, limit):
    result = 1
    for _ in range(exponent):
        if result > limit // base:
            return limit + 1
        result *= base
    return result

def sieve_primes(limit):
    if limit < 2: return []
    is_prime = bytearray([1]) * (limit + 1)
    is_prime[0] = is_prime[1] = 0
    
    for p in range(2, int(math.sqrt(limit)) + 1):
        if is_prime[p]:
            for m in range(p * p, limit + 1, p):
                is_prime[m] = 0
                
    return [p for p in range(2, limit + 1) if is_prime[p]]

def primes_1mod4_up_to(limit):
    return [p for p in sieve_primes(limit) if p % 4 == 1]

def build_spf(max_n):
    spf = bytearray(max_n + 1)
    # Use array of integers since max_n could be up to 10^7
    spf = [0] * (max_n + 1)
    primes = []
    
    for i in range(2, max_n + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        for p in primes:
            m = p * i
            if m > max_n:
                break
            spf[m] = p
            if p == spf[i]:
                break
    return spf

def tail_fits(ordered_exponents, next_pos, next_prime_idx, remaining_limit, primes_1mod4):
    idx = next_prime_idx
    for pos in range(next_pos, len(ordered_exponents)):
        if idx >= len(primes_1mod4):
            return False
        p = primes_1mod4[idx]
        pe = pow_with_limit(p, ordered_exponents[pos], remaining_limit)
        if pe > remaining_limit:
            return False
        remaining_limit //= pe
        idx += 1
    return True

def enumerate_core_for_order(ordered_exponents, pos, start_prime_idx, current_product, limit, primes_1mod4, out):
    if pos == len(ordered_exponents):
        out.append(current_product)
        return
        
    remaining_slots = len(ordered_exponents) - pos
    if start_prime_idx + remaining_slots > len(primes_1mod4):
        return
        
    head_limit = limit // current_product
    for i in range(start_prime_idx, len(primes_1mod4) - remaining_slots + 1):
        p = primes_1mod4[i]
        pe = pow_with_limit(p, ordered_exponents[pos], head_limit)
        if pe > head_limit:
            break
            
        next_product = current_product * pe
        if pos + 1 < len(ordered_exponents):
            if not tail_fits(ordered_exponents, pos + 1, i + 1, limit // next_product, primes_1mod4):
                break
                
        enumerate_core_for_order(ordered_exponents, pos + 1, i + 1, next_product, limit, primes_1mod4, out)

def enumerate_core_numbers(limit, primes_1mod4):
    exponent_patterns = [
        [52],
        [17, 1],
        [10, 2],
        [7, 3],
        [3, 2, 1],
    ]
    
    core_numbers = []
    for pattern in exponent_patterns:
        for perm in sorted(set(itertools.permutations(pattern))):
            enumerate_core_for_order(perm, 0, 0, 1, limit, primes_1mod4, core_numbers)
            
    return sorted(list(set(core_numbers)))

def build_allowed_multiplier_prefix(max_multiplier):
    prefix = [0] * (max_multiplier + 1)
    if max_multiplier == 0:
        return prefix
        
    spf = build_spf(max_multiplier)
    valid = bytearray(max_multiplier + 1)
    valid[1] = 1
    
    for n in range(2, max_multiplier + 1):
        p = spf[n]
        m = n // p
        allowed = (p == 2) or ((p % 4) == 3)
        valid[n] = 1 if (allowed and valid[m]) else 0
        
    running = 0
    for n in range(1, max_multiplier + 1):
        if valid[n]:
            running += n
        prefix[n] = running
        
    return prefix

def accumulate_range(args):
    core_numbers, begin, end, limit, allowed_prefix = args
    local_sum = 0
    for i in range(begin, end):
        core = core_numbers[i]
        max_mult = limit // core
        local_sum += core * allowed_prefix[max_mult]
    return local_sum

def solve_sum(limit, allow_multithreading=True):
    prime_bound = max(200, limit // 21125 + 100)
    primes_1mod4 = primes_1mod4_up_to(prime_bound)
    core_numbers = enumerate_core_numbers(limit, primes_1mod4)
    
    if not core_numbers:
        return 0
        
    min_core = core_numbers[0]
    max_multiplier = limit // min_core
    allowed_prefix = build_allowed_multiplier_prefix(max_multiplier)
    
    threads = multiprocessing.cpu_count() if allow_multithreading else 1
    threads = min(threads, max(1, len(core_numbers)))
    
    if threads == 1:
        return accumulate_range((core_numbers, 0, len(core_numbers), limit, allowed_prefix))
        
    pool_args = []
    for t in range(threads):
        begin = len(core_numbers) * t // threads
        end = len(core_numbers) * (t + 1) // threads
        pool_args.append((core_numbers, begin, end, limit, allowed_prefix))
        
    with multiprocessing.Pool(threads) as pool:
        results = pool.map(accumulate_range, pool_args)
        
    return sum(results)

def solve():
    return str(solve_sum(100000000000))

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler233 {
    static final long LIMIT = 100_000_000_000L;

    static long powWithLimit(long base, int exp, long limit) {
        long result = 1;
        for (int i = 0; i < exp; ++i) {
            if (result > limit / base) {
                return limit + 1;
            }
            result *= base;
        }
        return result;
    }

    static int[] primes1Mod4UpTo(int limit) {
        if (limit < 2)
            return new int[0];
        byte[] isPrime = new byte[limit + 1];
        Arrays.fill(isPrime, (byte) 1);
        isPrime[0] = isPrime[1] = 0;
        int root = (int) Math.sqrt(limit);
        for (int p = 2; p <= root; ++p) {
            if (isPrime[p] == 1) {
                for (long m = (long) p * p; m <= limit; m += p) {
                    isPrime[(int) m] = 0;
                }
            }
        }
        int count = 0;
        for (int p = 2; p <= limit; ++p) {
            if (isPrime[p] == 1 && p % 4 == 1)
                count++;
        }
        int[] primes = new int[count];
        int idx = 0;
        for (int p = 2; p <= limit; ++p) {
            if (isPrime[p] == 1 && p % 4 == 1)
                primes[idx++] = p;
        }
        return primes;
    }

    static int[] buildSpf(int maxN) {
        int[] spf = new int[maxN + 1];
        int[] primes = new int[maxN / 10 + 100];
        int numPrimes = 0;

        for (int i = 2; i <= maxN; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                if (numPrimes < primes.length) {
                    primes[numPrimes++] = i;
                }
            }
            for (int j = 0; j < numPrimes; ++j) {
                int p = primes[j];
                long m = (long) p * i;
                if (m > maxN)
                    break;
                spf[(int) m] = p;
                if (p == spf[i])
                    break;
            }
        }
        return spf;
    }

    static boolean tailFits(int[] orderedExponents, int nextPos, int nextPrimeIdx, long remainingLimit,
            int[] primes1Mod4) {
        int idx = nextPrimeIdx;
        for (int pos = nextPos; pos < orderedExponents.length; ++pos, ++idx) {
            if (idx >= primes1Mod4.length)
                return false;
            long p = primes1Mod4[idx];
            long pe = powWithLimit(p, orderedExponents[pos], remainingLimit);
            if (pe > remainingLimit)
                return false;
            remainingLimit /= pe;
        }
        return true;
    }

    static void enumerateCoreForOrder(int[] orderedExponents, int pos, int startPrimeIdx, long currentProduct,
            long limit, int[] primes1Mod4, List<Long> out) {
        if (pos == orderedExponents.length) {
            out.add(currentProduct);
            return;
        }
        int remainingSlots = orderedExponents.length - pos;
        if (startPrimeIdx + remainingSlots > primes1Mod4.length)
            return;

        long headLimit = limit / currentProduct;
        for (int i = startPrimeIdx; i + remainingSlots <= primes1Mod4.length; ++i) {
            long p = primes1Mod4[i];
            long pe = powWithLimit(p, orderedExponents[pos], headLimit);
            if (pe > headLimit)
                break;

            long nextProduct = currentProduct * pe;
            if (pos + 1 < orderedExponents.length) {
                if (!tailFits(orderedExponents, pos + 1, i + 1, limit / nextProduct, primes1Mod4)) {
                    break;
                }
            }
            enumerateCoreForOrder(orderedExponents, pos + 1, i + 1, nextProduct, limit, primes1Mod4, out);
        }
    }

    static void swap(int[] arr, int i, int j) {
        int t = arr[i];
        arr[i] = arr[j];
        arr[j] = t;
    }

    static void reverse(int[] arr, int i, int j) {
        while (i < j)
            swap(arr, i++, j--);
    }

    static boolean nextPermutation(int[] arr) {
        int i = arr.length - 2;
        while (i >= 0 && arr[i] >= arr[i + 1])
            i--;
        if (i < 0)
            return false;
        int j = arr.length - 1;
        while (arr[j] <= arr[i])
            j--;
        swap(arr, i, j);
        reverse(arr, i + 1, arr.length - 1);
        return true;
    }

    static long[] enumerateCoreNumbers(long limit, int[] primes1Mod4) {
        int[][] exponentPatterns = {
                { 52 }, { 1, 17 }, { 2, 10 }, { 3, 7 }, { 1, 2, 3 }
        };

        List<Long> coreList = new ArrayList<>();
        for (int[] pattern : exponentPatterns) {
            Arrays.sort(pattern);
            do {
                enumerateCoreForOrder(pattern, 0, 0, 1L, limit, primes1Mod4, coreList);
            } while (nextPermutation(pattern));
        }

        Collections.sort(coreList);
        int uniqueCount = 0;
        for (int i = 0; i < coreList.size(); ++i) {
            if (i == 0 || !coreList.get(i).equals(coreList.get(i - 1))) {
                uniqueCount++;
            }
        }
        long[] coreNumbers = new long[uniqueCount];
        int idx = 0;
        for (int i = 0; i < coreList.size(); ++i) {
            if (i == 0 || !coreList.get(i).equals(coreList.get(i - 1))) {
                coreNumbers[idx++] = coreList.get(i);
            }
        }
        return coreNumbers;
    }

    static long[] buildAllowedMultiplierPrefix(int maxMultiplier) {
        long[] prefix = new long[maxMultiplier + 1];
        if (maxMultiplier == 0)
            return prefix;

        int[] spf = buildSpf(maxMultiplier);
        byte[] valid = new byte[maxMultiplier + 1];
        valid[1] = 1;

        for (int n = 2; n <= maxMultiplier; ++n) {
            int p = spf[n];
            int m = n / p;
            boolean allowed = (p == 2) || ((p % 4) == 3);
            if (allowed && valid[m] == 1) {
                valid[n] = 1;
            }
        }

        long running = 0;
        for (int n = 1; n <= maxMultiplier; ++n) {
            if (valid[n] == 1) {
                running += n;
            }
            prefix[n] = running;
        }
        return prefix;
    }

    static long accumulateRange(long[] coreNumbers, int begin, int end, long limit, long[] allowedPrefix) {
        long local = 0;
        for (int i = begin; i < end; ++i) {
            long core = coreNumbers[i];
            int maxMult = (int) (limit / core);
            local += core * allowedPrefix[maxMult];
        }
        return local;
    }

    public static String solve() {
        int primeBound = (int) Math.max(200L, LIMIT / 21125 + 100);
        int[] primes1Mod4 = primes1Mod4UpTo(primeBound);
        long[] coreNumbers = enumerateCoreNumbers(LIMIT, primes1Mod4);

        if (coreNumbers.length == 0)
            return "0";

        long minCore = coreNumbers[0];
        int maxMultiplier = (int) (LIMIT / minCore);
        long[] allowedPrefix = buildAllowedMultiplierPrefix(maxMultiplier);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        threads = Math.min(threads, coreNumbers.length);

        if (threads <= 1) {
            return String.valueOf(accumulateRange(coreNumbers, 0, coreNumbers.length, LIMIT, allowedPrefix));
        }

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int t = 0; t < threads; ++t) {
            final int begin = coreNumbers.length * t / threads;
            final int end = coreNumbers.length * (t + 1) / threads;
            futures.add(executor.submit(() -> accumulateRange(coreNumbers, begin, end, LIMIT, allowedPrefix)));
        }

        long total = 0;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
            }
        }
        executor.shutdown();

        return String.valueOf(total);
    }

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