Problem 229: Four Representations Using Squares

Project Euler Problem 229 Solution

EulerSolve provides an optimized solution for Project Euler Problem 229, Four Representations Using Squares, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count integers $n \le 2\times 10^9$ that can all be written in the four forms $$n=a^2+b^2=a^2+2b^2=a^2+3b^2=a^2+7b^2,$$ with positive integers $a$ and $b$ in each representation. The implementations do not search directly over all pairs $(a,b)$. Instead they first count a slightly larger set in which zero is temporarily allowed, classify those integers by their squarefree kernel, and then repair the overcount caused by perfect squares. Mathematical Approach Let $L=2\times 10^9$ and let $D=\{1,2,3,7\}$. For the main counting argument it is convenient to relax the condition and first consider integers $n$ such that for every $d\in D$ there exist $a_d,b_d\ge 0$ with $$n=a_d^2+d\,b_d^2.$$ The strict positivity condition will be restored afterward. Squarefree Kernels Control the Relaxed Problem Every positive integer has a unique decomposition $$n=q\,m^2,$$ where $q$ is squarefree. The implementations take as their number-theoretic starting point the characterization that simultaneous representability in all four relaxed forms depends only on this squarefree kernel....

Detailed mathematical approach

Problem Summary

We must count integers $n \le 2\times 10^9$ that can all be written in the four forms

$$n=a^2+b^2=a^2+2b^2=a^2+3b^2=a^2+7b^2,$$

with positive integers $a$ and $b$ in each representation.

The implementations do not search directly over all pairs $(a,b)$. Instead they first count a slightly larger set in which zero is temporarily allowed, classify those integers by their squarefree kernel, and then repair the overcount caused by perfect squares.

Mathematical Approach

Let $L=2\times 10^9$ and let $D=\{1,2,3,7\}$. For the main counting argument it is convenient to relax the condition and first consider integers $n$ such that for every $d\in D$ there exist $a_d,b_d\ge 0$ with

$$n=a_d^2+d\,b_d^2.$$

The strict positivity condition will be restored afterward.

Squarefree Kernels Control the Relaxed Problem

Every positive integer has a unique decomposition

$$n=q\,m^2,$$

where $q$ is squarefree. The implementations take as their number-theoretic starting point the characterization that simultaneous representability in all four relaxed forms depends only on this squarefree kernel. More precisely, the relaxed condition holds exactly when every prime factor of $q$ lies in the residue classes

$$p\equiv 1,\ 25,\ 121 \pmod{168}.$$

So the admissible kernels are precisely

$$q=\prod_{i=1}^{r} p_i,$$

where the $p_i$ are distinct primes from those three classes; $q=1$ is the empty product. Once such a kernel is fixed, every number of the form $q\,m^2$ belongs to the relaxed set, and once the kernel contains any other prime, the number can be discarded immediately.

Counting Admissible Kernels

For a fixed admissible kernel $q$, the valid integers are exactly the numbers $q\,m^2\le L$. Their count is therefore

$$\left\lfloor \sqrt{\frac{L}{q}} \right\rfloor.$$

If $R(L)$ denotes the relaxed count, then

$$R(L)=\sum_{\substack{q\le L\\ q\text{ admissible}}}\left\lfloor \sqrt{\frac{L}{q}} \right\rfloor.$$

The implementations split this sum into three disjoint cases.

The kernel $q=1$ contributes $\lfloor\sqrt{L}\rfloor$. A singleton kernel $q=p$ contributes $\lfloor\sqrt{L/p}\rfloor$ for each good prime $p$. All remaining kernels contain at least two distinct good primes, so the implementations enumerate their products by depth-first search and add the same floor term for each feasible product.

An important pruning rule comes from the smallest good prime. Since the first prime with residue $1$, $25$, or $121$ modulo $168$ is $193$, any admissible kernel with at least two prime factors is at least $193^2$. That is why the multi-prime stage can stop very early for small partial products.

Why Perfect Squares Need a Separate Correction

The relaxed count is not yet the answer, because the original problem requires positive $a$ and $b$, not merely nonnegative ones. The only systematic overcount comes from perfect squares.

Indeed, every square $m^2$ is automatically included in the relaxed set via the trivial identities

$$m^2=m^2+d\cdot 0^2 \qquad (d\in D).$$

For a nonsquare, the degenerate case $b=0$ cannot occur. And if $a=0$ for $d=2$, $3$, or $7$, then the squarefree kernel would contain $2$, $3$, or $7$, which is impossible for an admissible kernel. So the positivity repair can be written as

$$A(L)=R(L)-\lfloor\sqrt{L}\rfloor+G(L),$$

where $G(L)$ counts the squares $m^2\le L$ that really do admit positive representations for all four coefficients.

Turning Square Representations into a Divisor Problem

To compute $G(L)$, fix $d\in D$ and ask when a square satisfies

$$m^2=a^2+d\,b^2$$

with $a,b>0$. Rearranging gives

$$(m-a)(m+a)=d\,b^2.$$

If we set

$$u=m-a,\qquad v=m+a,$$

then we obtain

$$uv=d\,b^2,\qquad u<v,\qquad u\equiv v\pmod 2,$$

and conversely any factor pair $(u,v)$ satisfying these conditions yields

$$m=\frac{u+v}{2},\qquad a=\frac{v-u}{2}.$$

That turns the search for positive square representations into a divisor enumeration problem. For each fixed $d$ and each positive $b$, factor $d\,b^2$, generate all divisors $u$, define $v=(d\,b^2)/u$, keep only pairs with the same parity and $v>u$, and mark the resulting $m$. Doing this separately for $d=1,2,3,7$ produces four boolean tables over $m\le \sqrt{L}$, and their intersection is exactly the set counted by $G(L)$.

Worked Example: The Square $3600$

The number $3600=60^2$ illustrates why the correction term is necessary. Its squarefree kernel is $q=1$, so it is counted automatically in the relaxed total $R(L)$.

It also survives the positivity correction, because it genuinely has positive representations in all four forms:

$$3600=48^2+36^2=20^2+2\cdot 40^2=30^2+3\cdot 30^2=45^2+7\cdot 15^2.$$

So $3600$ is first counted among all squares through the kernel $q=1$, then removed by the blanket subtraction of $\lfloor\sqrt{L}\rfloor$, and finally added back because it belongs to $G(L)$.

How the Code Works

The C++, Python, and Java implementations follow exactly the split described above. First they compute the relaxed quantity $R(L)$. The singleton-kernel contribution is obtained by scanning primes up to $L$ with a segmented sieve and keeping only primes in the three admissible residue classes modulo $168$. Each such prime $p$ contributes $\left\lfloor \sqrt{L/p} \right\rfloor$.

The multi-prime part uses the ordered list of good primes up to $L/193$. A depth-first search grows products of distinct good primes in increasing order, stops as soon as the next multiplication would exceed $L$, and adds $\left\lfloor \sqrt{L/q} \right\rfloor$ whenever the current product $q$ already contains at least two primes. In this way every admissible multi-prime kernel is visited exactly once.

The square correction then builds a smallest-prime-factor table up to $\sqrt{L}$. For each $d\in\{1,2,3,7\}$ and each positive $b\le \sqrt{L/d}$, the implementation factors $d\,b^2$, generates all divisors, reconstructs every valid $m$ from the factor pairs $(u,v)$, and marks it. The final answer is the relaxed kernel count minus all squares plus the number of $m$ marked in every one of the four tables.

The language-specific execution strategy differs slightly, but the mathematics is the same. The C++ and Java implementations parallelize the segmented prime sweep across workers, while the Python implementation uses multiple processes for that same stage.

Complexity Analysis

The dominant large-scale cost is the singleton-kernel stage, which performs a segmented prime sweep over the interval up to $L$ while sieving with the base primes up to $\sqrt{L}$. In asymptotic terms this is standard segmented-sieve work over a range of length $L$, while memory stays bounded by the segment size and the small base-prime list.

The multi-prime DFS is much smaller in practice, because admissible kernels grow quickly once distinct good primes are multiplied together, and the recursion prunes immediately when the running product would exceed $L$. The square-correction phase only works up to $\sqrt{L}$, using a smallest-prime-factor table and divisor generation for numbers of the form $d\,b^2$.

Memory usage is modest throughout: one segment buffer for the sieve, prime-factor data up to $\sqrt{L}$, and four mark arrays of length $\lfloor\sqrt{L}\rfloor$. The whole algorithm is efficient because it replaces a huge search in the $(a,b)$-plane by kernel counting plus a focused correction on squares.

Footnotes and References

Problem page: https://projecteuler.net/problem=229
Binary quadratic form: Wikipedia - Binary quadratic form
Sum of two squares theorem: Wikipedia - Sum of two squares theorem
Squarefree integer: Wikipedia - Squarefree integer
Segmented sieve: cp-algorithms - Sieve of Eratosthenes

Problem 229 source code

C++

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

namespace {

using u64 = std::uint64_t;

constexpr u64 kDefaultLimit = 2000000000ULL;
constexpr u64 kMinGoodPrime = 193ULL;

struct Options {
    u64 limit = kDefaultLimit;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    bool run_checkpoints = true;
};

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

bool parse_u64_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

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

    value = parsed;
    return true;
}

bool parse_int_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    long long parsed = 0;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10LL + static_cast<long long>(ch - '0');
        if (parsed > static_cast<long long>(std::numeric_limits<int>::max())) {
            return false;
        }
    }

    value = static_cast<int>(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 (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }

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

    if (options.threads <= 0) {
        options.threads = 1;
    }

    return true;
}

std::vector<int> sieve_primes_up_to(const int n) {
    if (n < 2) {
        return {};
    }

    std::vector<char> is_prime(static_cast<std::size_t>(n) + 1U, 1);
    is_prime[0] = 0;
    is_prime[1] = 0;

    const int root = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    for (int p = 2; p <= root; ++p) {
        if (!is_prime[p]) {
            continue;
        }
        const int start = p * p;
        for (int x = start; x <= n; x += p) {
            is_prime[x] = 0;
        }
    }

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

bool is_good_prime(const u64 p) {
    if (p <= 7ULL) {
        return false;
    }
    const u64 residue = p % 168ULL;
    return residue == 1ULL || residue == 25ULL || residue == 121ULL;
}

u64 singleton_sum_in_range(const u64 low,
                           const u64 high,
                           const u64 limit,
                           const std::vector<int>& base_primes) {
    if (low > high || high < 2ULL) {
        return 0ULL;
    }

    constexpr u64 segment_size = 1ULL << 20;
    u64 total = 0ULL;

    for (u64 segment_low = low; segment_low <= high;) {
        const u64 segment_high =
            std::min(segment_low + segment_size - 1ULL, high);
        const std::size_t length =
            static_cast<std::size_t>(segment_high - segment_low + 1ULL);

        std::vector<char> is_prime(length, 1);

        for (const int p_int : base_primes) {
            const u64 p = static_cast<u64>(p_int);
            const u64 p2 = p * p;
            if (p2 > segment_high) {
                break;
            }

            u64 start = (segment_low + p - 1ULL) / p;
            start *= p;
            if (start < p2) {
                start = p2;
            }

            for (u64 x = start; x <= segment_high; x += p) {
                is_prime[static_cast<std::size_t>(x - segment_low)] = 0;
            }
        }

        if (segment_low == 0ULL) {
            is_prime[0] = 0;
            if (length > 1U) {
                is_prime[1] = 0;
            }
        } else if (segment_low == 1ULL) {
            is_prime[0] = 0;
        }

        for (u64 value = std::max<u64>(2ULL, segment_low); value <= segment_high;
             ++value) {
            if (!is_prime[static_cast<std::size_t>(value - segment_low)]) {
                continue;
            }
            if (!is_good_prime(value)) {
                continue;
            }
            total += isqrt_u64(limit / value);
        }

        if (segment_high == std::numeric_limits<u64>::max()) {
            break;
        }
        segment_low = segment_high + 1ULL;
    }

    return total;
}

u64 sum_singleton_kernels(const u64 limit, const int threads) {
    if (limit < 11ULL) {
        return 0ULL;
    }

    const int base_bound = static_cast<int>(isqrt_u64(limit));
    const std::vector<int> base_primes = sieve_primes_up_to(base_bound);

    const u64 first = 2ULL;
    const u64 last = limit;
    const u64 count = last - first + 1ULL;

    int worker_count = std::max(1, threads);
    if (static_cast<u64>(worker_count) > count) {
        worker_count = static_cast<int>(count);
    }

    std::vector<u64> partial(static_cast<std::size_t>(worker_count), 0ULL);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(worker_count));

    const u64 block = (count + static_cast<u64>(worker_count) - 1ULL) /
                      static_cast<u64>(worker_count);

    for (int tid = 0; tid < worker_count; ++tid) {
        const u64 low = first + static_cast<u64>(tid) * block;
        if (low > last) {
            break;
        }
        const u64 high = std::min(last, low + block - 1ULL);

        workers.emplace_back([&, tid, low, high]() {
            partial[static_cast<std::size_t>(tid)] =
                singleton_sum_in_range(low, high, limit, base_primes);
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

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

std::vector<int> collect_good_primes_up_to(const int n) {
    if (n < 11) {
        return {};
    }

    const std::vector<int> primes = sieve_primes_up_to(n);
    std::vector<int> good;
    good.reserve(primes.size() / 12U);

    for (const int p : primes) {
        if (is_good_prime(static_cast<u64>(p))) {
            good.push_back(p);
        }
    }

    return good;
}

u64 sum_multi_prime_kernels_dfs(const std::vector<int>& good_primes,
                                const u64 limit,
                                const std::size_t start_index,
                                const u64 current_product,
                                const int picked_count) {
    u64 total = 0ULL;

    for (std::size_t i = start_index; i < good_primes.size(); ++i) {
        const u64 p = static_cast<u64>(good_primes[i]);
        if (current_product > limit / p) {
            break;
        }

        const u64 next_product = current_product * p;
        const int next_picked = picked_count + 1;

        if (next_picked >= 2) {
            total += isqrt_u64(limit / next_product);
        }

        total += sum_multi_prime_kernels_dfs(good_primes,
                                             limit,
                                             i + 1,
                                             next_product,
                                             next_picked);
    }

    return total;
}

u64 sum_multi_prime_kernels(const u64 limit) {
    if (limit < kMinGoodPrime * kMinGoodPrime) {
        return 0ULL;
    }

    const u64 max_factor = limit / kMinGoodPrime;
    if (max_factor > static_cast<u64>(std::numeric_limits<int>::max())) {
        throw std::runtime_error("max_factor does not fit into int");
    }

    const std::vector<int> good_primes =
        collect_good_primes_up_to(static_cast<int>(max_factor));

    return sum_multi_prime_kernels_dfs(good_primes, limit, 0U, 1ULL, 0);
}

u64 count_integer_representable(const u64 limit, const int threads) {
    const u64 square_kernel_contribution = isqrt_u64(limit);
    const u64 singleton_contribution = sum_singleton_kernels(limit, threads);
    const u64 multi_prime_contribution = sum_multi_prime_kernels(limit);

    return square_kernel_contribution + singleton_contribution +
           multi_prime_contribution;
}

std::vector<int> build_spf_table(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n) + 1U, 0);
    for (int i = 0; i <= n; ++i) {
        spf[static_cast<std::size_t>(i)] = i;
    }

    if (n >= 0) {
        spf[0] = 0;
    }
    if (n >= 1) {
        spf[1] = 1;
    }

    const int root = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    for (int i = 2; i <= root; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (int x = i * i; x <= n; x += i) {
            if (spf[static_cast<std::size_t>(x)] == x) {
                spf[static_cast<std::size_t>(x)] = i;
            }
        }
    }

    return spf;
}

void add_factor(std::vector<std::pair<int, int>>& factors,
                const int prime,
                const int exponent) {
    for (auto& [p, e] : factors) {
        if (p == prime) {
            e += exponent;
            return;
        }
    }
    factors.push_back({prime, exponent});
}

std::vector<std::pair<int, int>> factorize_dy2(const int y,
                                                const int d,
                                                const std::vector<int>& spf) {
    std::vector<std::pair<int, int>> factors;

    int value = y;
    while (value > 1) {
        const int p = spf[static_cast<std::size_t>(value)];
        int exponent = 0;
        while (value % p == 0) {
            value /= p;
            ++exponent;
        }
        add_factor(factors, p, 2 * exponent);
    }

    int remaining_d = d;
    for (int p = 2; static_cast<long long>(p) * p <= remaining_d; ++p) {
        int exponent = 0;
        while (remaining_d % p == 0) {
            remaining_d /= p;
            ++exponent;
        }
        if (exponent > 0) {
            add_factor(factors, p, exponent);
        }
    }
    if (remaining_d > 1) {
        add_factor(factors, remaining_d, 1);
    }

    return factors;
}

void build_divisors_recursive(const std::vector<std::pair<int, int>>& factors,
                              const std::size_t index,
                              const u64 current,
                              std::vector<u64>& out) {
    if (index == factors.size()) {
        out.push_back(current);
        return;
    }

    const int prime = factors[index].first;
    const int exponent = factors[index].second;

    u64 value = current;
    for (int i = 0; i <= exponent; ++i) {
        build_divisors_recursive(factors, index + 1U, value, out);
        value *= static_cast<u64>(prime);
    }
}

std::vector<char> square_roots_with_positive_rep(const u64 limit,
                                                 const int d,
                                                 const std::vector<int>& spf) {
    const int max_root = static_cast<int>(isqrt_u64(limit));
    const int max_y = static_cast<int>(isqrt_u64(limit / static_cast<u64>(d)));

    std::vector<char> mark(static_cast<std::size_t>(max_root) + 1U, 0);

    std::vector<u64> divisors;
    divisors.reserve(512);

    for (int y = 1; y <= max_y; ++y) {
        const u64 D = static_cast<u64>(d) * static_cast<u64>(y) *
                      static_cast<u64>(y);

        const std::vector<std::pair<int, int>> factors = factorize_dy2(y, d, spf);

        divisors.clear();
        build_divisors_recursive(factors, 0U, 1ULL, divisors);

        for (const u64 u : divisors) {
            const u64 v = D / u;
            if (u > v) {
                continue;
            }
            if (((u + v) & 1ULL) != 0ULL) {
                continue;
            }
            if (v <= u) {
                continue;
            }

            const u64 m = (u + v) / 2ULL;
            if (m <= static_cast<u64>(max_root)) {
                mark[static_cast<std::size_t>(m)] = 1;
            }
        }
    }

    return mark;
}

u64 count_good_squares(const u64 limit) {
    const int max_root = static_cast<int>(isqrt_u64(limit));
    const std::vector<int> spf = build_spf_table(max_root);

    const std::vector<char> d1 = square_roots_with_positive_rep(limit, 1, spf);
    const std::vector<char> d2 = square_roots_with_positive_rep(limit, 2, spf);
    const std::vector<char> d3 = square_roots_with_positive_rep(limit, 3, spf);
    const std::vector<char> d7 = square_roots_with_positive_rep(limit, 7, spf);

    u64 count = 0ULL;
    for (int m = 1; m <= max_root; ++m) {
        if (d1[static_cast<std::size_t>(m)] && d2[static_cast<std::size_t>(m)] &&
            d3[static_cast<std::size_t>(m)] && d7[static_cast<std::size_t>(m)]) {
            ++count;
        }
    }
    return count;
}

u64 count_special_numbers(const u64 limit, const int threads) {
    const u64 integer_representable = count_integer_representable(limit, threads);
    const u64 total_squares = isqrt_u64(limit);
    const u64 good_squares = count_good_squares(limit);

    return integer_representable - total_squares + good_squares;
}

bool has_positive_rep_bruteforce(const u64 n, const int d) {
    for (u64 b = 1ULL;; ++b) {
        const u64 db2 = static_cast<u64>(d) * b * b;
        if (db2 >= n) {
            break;
        }

        const u64 rem = n - db2;
        const u64 a = isqrt_u64(rem);
        if (a > 0ULL && a * a == rem) {
            return true;
        }
    }
    return false;
}

u64 brute_force_count(const int limit) {
    const std::array<int, 4> ds = {1, 2, 3, 7};
    std::array<std::vector<char>, 4> mark;

    for (int i = 0; i < 4; ++i) {
        mark[i].assign(static_cast<std::size_t>(limit) + 1U, 0);
        const int d = ds[i];

        for (int b = 1;; ++b) {
            const u64 db2 = static_cast<u64>(d) * static_cast<u64>(b) *
                            static_cast<u64>(b);
            if (db2 > static_cast<u64>(limit)) {
                break;
            }
            const int max_a =
                static_cast<int>(isqrt_u64(static_cast<u64>(limit) - db2));
            for (int a = 1; a <= max_a; ++a) {
                const int n = static_cast<int>(static_cast<u64>(a) *
                                                   static_cast<u64>(a) +
                                               db2);
                mark[i][static_cast<std::size_t>(n)] = 1;
            }
        }
    }

    u64 count = 0ULL;
    for (int n = 1; n <= limit; ++n) {
        if (mark[0][static_cast<std::size_t>(n)] &&
            mark[1][static_cast<std::size_t>(n)] &&
            mark[2][static_cast<std::size_t>(n)] &&
            mark[3][static_cast<std::size_t>(n)]) {
            ++count;
        }
    }

    return count;
}

void run_checkpoints(const int threads) {
    const std::array<u64, 2> examples = {3600ULL, 88201ULL};
    for (const u64 n : examples) {
        for (const int d : {1, 2, 3, 7}) {
            if (!has_positive_rep_bruteforce(n, d)) {
                throw std::runtime_error("Example validation failed");
            }
        }
    }

    constexpr int small_limit = 200000;
    const u64 fast_small =
        count_special_numbers(static_cast<u64>(small_limit), threads);
    const u64 brute_small = brute_force_count(small_limit);
    if (fast_small != brute_small) {
        throw std::runtime_error("Small brute-force checkpoint failed");
    }

    constexpr u64 stated_limit = 10000000ULL;
    constexpr u64 stated_value = 75373ULL;
    const u64 fast_stated = count_special_numbers(stated_limit, threads);
    if (fast_stated != stated_value) {
        throw std::runtime_error("10^7 checkpoint failed");
    }
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        std::cerr << "Usage: ./Euler229 [--limit=N] [--threads=T] [--skip-checkpoints]\n";
        return 1;
    }

    try {
        if (options.run_checkpoints) {
            run_checkpoints(options.threads);
        }

        const u64 answer = count_special_numbers(options.limit, options.threads);
        std::cout << answer << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

from __future__ import annotations

import bisect
import math
import multiprocessing as mp
import os
from concurrent.futures import ProcessPoolExecutor


G_LIMIT = 0
G_BASE_PRIMES: list[int] = []


def isqrt(n: int) -> int:
    return math.isqrt(n)


def sieve_primes_up_to(n: int) -> list[int]:
    if n < 2:
        return []
    is_p = bytearray(b"\x01" * (n + 1))
    is_p[0] = 0
    is_p[1] = 0
    root = isqrt(n)
    for p in range(2, root + 1):
        if is_p[p]:
            start = p * p
            is_p[start::p] = bytearray(((n - start) // p) + 1)
    return [p for p in range(2, n + 1) if is_p[p]]


def is_good_prime(p: int) -> bool:
    if p <= 7:
        return False
    residue = p % 168
    return residue == 1 or residue == 25 or residue == 121


def singleton_sum_in_range(low: int, high: int, limit: int, base_primes: list[int]) -> int:
    if low > high or high < 2:
        return 0

    segment_size = 1 << 20
    total = 0
    seg_low = low

    while seg_low <= high:
        seg_high = min(seg_low + segment_size - 1, high)
        length = seg_high - seg_low + 1
        is_prime = bytearray(b"\x01" * length)

        for p in base_primes:
            p2 = p * p
            if p2 > seg_high:
                break
            start = ((seg_low + p - 1) // p) * p
            if start < p2:
                start = p2
            is_prime[start - seg_low : length : p] = bytearray(((seg_high - start) // p) + 1)

        if seg_low <= 1:
            bound = min(2, length)
            for i in range(bound):
                if seg_low + i <= 1:
                    is_prime[i] = 0

        for idx, flag in enumerate(is_prime):
            if not flag:
                continue
            value = seg_low + idx
            if is_good_prime(value):
                total += isqrt(limit // value)

        seg_low = seg_high + 1

    return total


def init_singleton_worker(limit: int, base_primes: list[int]) -> None:
    global G_LIMIT, G_BASE_PRIMES
    G_LIMIT = limit
    G_BASE_PRIMES = base_primes


def singleton_worker(bounds: tuple[int, int]) -> int:
    low, high = bounds
    return singleton_sum_in_range(low, high, G_LIMIT, G_BASE_PRIMES)


def sum_singleton_kernels(limit: int, base_primes: list[int], threads: int) -> int:
    first = 2
    last = limit
    count = last - first + 1
    workers = max(1, min(threads, count))
    if workers == 1:
        return singleton_sum_in_range(first, last, limit, base_primes)

    block = (count + workers - 1) // workers
    jobs = []
    for tid in range(workers):
        low = first + tid * block
        if low > last:
            break
        high = min(last, low + block - 1)
        jobs.append((low, high))

    ctx = mp.get_context("fork") if "fork" in mp.get_all_start_methods() else mp.get_context()
    total = 0
    with ProcessPoolExecutor(
        max_workers=workers,
        mp_context=ctx,
        initializer=init_singleton_worker,
        initargs=(limit, base_primes),
    ) as pool:
        for partial in pool.map(singleton_worker, jobs):
            total += partial
    return total


def build_spf_table(n: int) -> list[int]:
    spf = list(range(n + 1))
    root = isqrt(n)
    for i in range(2, root + 1):
        if spf[i] == i:
            for x in range(i * i, n + 1, i):
                if spf[x] == x:
                    spf[x] = i
    return spf


def factorize_dy2(y: int, d: int, spf: list[int]) -> list[tuple[int, int]]:
    factors: dict[int, int] = {}
    val = y
    while val > 1:
        p = spf[val]
        count = 0
        while val % p == 0:
            val //= p
            count += 1
        factors[p] = factors.get(p, 0) + 2 * count

    rem_d = d
    p = 2
    while p * p <= rem_d:
        if rem_d % p == 0:
            count = 0
            while rem_d % p == 0:
                rem_d //= p
                count += 1
            factors[p] = factors.get(p, 0) + count
        p += 1
    if rem_d > 1:
        factors[rem_d] = factors.get(rem_d, 0) + 1

    return sorted(factors.items())


def build_divisors(factors: list[tuple[int, int]]) -> list[int]:
    divisors = [1]
    for p, e in factors:
        base = list(divisors)
        mul = 1
        for _ in range(e):
            mul *= p
            for d in base:
                divisors.append(d * mul)
    return divisors


def square_roots_with_positive_rep(limit: int, d: int, spf: list[int]) -> bytearray:
    max_root = isqrt(limit)
    max_y = isqrt(limit // d)
    mark = bytearray(max_root + 1)

    for y in range(1, max_y + 1):
        D = d * y * y
        factors = factorize_dy2(y, d, spf)
        divisors = build_divisors(factors)

        for u in divisors:
            v = D // u
            if u > v:
                continue
            if ((u + v) & 1) != 0:
                continue
            if v <= u:
                continue
            m = (u + v) // 2
            if m <= max_root:
                mark[m] = 1

    return mark


def solve() -> str:
    limit = 2_000_000_000
    threads = max(1, min(8, os.cpu_count() or 1))

    square_kernel = isqrt(limit)
    base_bound = isqrt(limit)
    base_primes = sieve_primes_up_to(base_bound)
    singleton = sum_singleton_kernels(limit, base_primes, threads)

    min_good = 193
    max_factor = limit // min_good if min_good <= limit else 0
    good_primes = [p for p in sieve_primes_up_to(max_factor) if is_good_prime(p)] if max_factor >= 11 else []

    def multi_dfs(start_index: int, current_product: int, picked: int) -> int:
        total = 0
        for i in range(start_index, len(good_primes)):
            p = good_primes[i]
            if current_product > limit // p:
                break
            next_product = current_product * p
            if picked + 1 >= 2:
                total += isqrt(limit // next_product)
            total += multi_dfs(i + 1, next_product, picked + 1)
        return total

    multi = multi_dfs(0, 1, 0)
    integer_representable = square_kernel + singleton + multi

    max_root = isqrt(limit)
    spf = build_spf_table(max_root)
    d1 = square_roots_with_positive_rep(limit, 1, spf)
    d2 = square_roots_with_positive_rep(limit, 2, spf)
    d3 = square_roots_with_positive_rep(limit, 3, spf)
    d7 = square_roots_with_positive_rep(limit, 7, spf)

    good_squares = sum(1 for m in range(1, max_root + 1) if d1[m] and d2[m] and d3[m] and d7[m])
    total_squares = isqrt(limit)
    answer = integer_representable - total_squares + good_squares
    return str(answer)


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

Java

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

public class Euler229 {

    static long isqrt(long n) {
        if (n == 0)
            return 0;
        long x = (long) Math.sqrt(n);
        while ((x + 1) * (x + 1) <= n)
            x++;
        while (x * x > n)
            x--;
        return x;
    }

    static int[] getBasePrimes(int n) {
        if (n < 2)
            return new int[0];
        byte[] isPrime = new byte[n + 1];
        Arrays.fill(isPrime, (byte) 1);
        isPrime[0] = isPrime[1] = 0;
        int root = (int) isqrt(n);
        for (int p = 2; p <= root; ++p) {
            if (isPrime[p] == 1) {
                for (int x = p * p; x <= n; x += p)
                    isPrime[x] = 0;
            }
        }
        int count = 0;
        for (int p = 2; p <= n; ++p)
            if (isPrime[p] == 1)
                count++;
        int[] primes = new int[count];
        int idx = 0;
        for (int p = 2; p <= n; ++p)
            if (isPrime[p] == 1)
                primes[idx++] = p;
        return primes;
    }

    static boolean isGoodPrime(long p) {
        if (p <= 7)
            return false;
        long residue = p % 168;
        return residue == 1 || residue == 25 || residue == 121;
    }

    static long singletonSumInRange(long low, long high, long limit, int[] basePrimes) {
        if (low > high || high < 2)
            return 0;
        int segmentSize = 1 << 20;
        long total = 0;

        for (long segmentLow = low; segmentLow <= high;) {
            long segmentHigh = Math.min(segmentLow + segmentSize - 1, high);
            int length = (int) (segmentHigh - segmentLow + 1);
            byte[] isPrime = new byte[length];
            Arrays.fill(isPrime, (byte) 1);

            for (int p : basePrimes) {
                long p2 = (long) p * p;
                if (p2 > segmentHigh)
                    break;
                long start = (segmentLow + p - 1) / p * p;
                if (start < p2)
                    start = p2;
                for (long x = start; x <= segmentHigh; x += p) {
                    isPrime[(int) (x - segmentLow)] = 0;
                }
            }

            if (segmentLow == 0) {
                isPrime[0] = 0;
                if (length > 1)
                    isPrime[1] = 0;
            } else if (segmentLow == 1) {
                isPrime[0] = 0;
            }

            for (long val = Math.max(2L, segmentLow); val <= segmentHigh; ++val) {
                if (isPrime[(int) (val - segmentLow)] == 1 && isGoodPrime(val)) {
                    total += isqrt(limit / val);
                }
            }
            segmentLow = segmentHigh + 1;
        }
        return total;
    }

    static long sumSingletonKernels(long limit, int threads) {
        if (limit < 11)
            return 0;
        int baseBound = (int) isqrt(limit);
        int[] basePrimes = getBasePrimes(baseBound);

        long last = limit;
        long count = last - 1;
        int workerCount = Math.min(threads, (int) Math.min(count, Integer.MAX_VALUE));
        long block = (count + workerCount - 1) / workerCount;

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

        for (int tid = 0; tid < workerCount; ++tid) {
            final long low = 2 + tid * block;
            if (low > last)
                break;
            final long high = Math.min(last, low + block - 1);
            futures.add(executor.submit(() -> singletonSumInRange(low, high, limit, basePrimes)));
        }

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

    static long sumMultiPrimeKernelsDfs(int[] goodPrimes, long limit, int startIndex, long currentProduct,
            int pickedCount) {
        long total = 0;
        for (int i = startIndex; i < goodPrimes.length; ++i) {
            long p = goodPrimes[i];
            if (currentProduct > limit / p)
                break;
            long nextProduct = currentProduct * p;
            int nextPicked = pickedCount + 1;
            if (nextPicked >= 2) {
                total += isqrt(limit / nextProduct);
            }
            total += sumMultiPrimeKernelsDfs(goodPrimes, limit, i + 1, nextProduct, nextPicked);
        }
        return total;
    }

    static long sumMultiPrimeKernels(long limit) {
        if (limit < 193L * 193L)
            return 0;
        int maxFactor = (int) (limit / 193);
        int[] basePrimes = getBasePrimes(maxFactor);

        int cnt = 0;
        for (int p : basePrimes)
            if (isGoodPrime(p))
                cnt++;
        int[] goodPrimes = new int[cnt];
        int idx = 0;
        for (int p : basePrimes)
            if (isGoodPrime(p))
                goodPrimes[idx++] = p;

        return sumMultiPrimeKernelsDfs(goodPrimes, limit, 0, 1L, 0);
    }

    static int[] buildSpfTable(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; ++i)
            spf[i] = i;
        int root = (int) isqrt(n);
        for (int i = 2; i <= root; ++i) {
            if (spf[i] == i) {
                for (int x = i * i; x <= n; x += i) {
                    if (spf[x] == x)
                        spf[x] = i;
                }
            }
        }
        return spf;
    }

    static class Factor {
        int p, e;

        Factor(int p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static List<Factor> factorizeDy2(int y, int d, int[] spf) {
        List<Factor> factors = new ArrayList<>();
        int val = y;
        while (val > 1) {
            int p = spf[val];
            int count = 0;
            while (val % p == 0) {
                val /= p;
                count++;
            }
            factors.add(new Factor(p, 2 * count));
        }

        int remD = d;
        for (int p = 2; p * p <= remD; ++p) {
            if (remD % p == 0) {
                int count = 0;
                while (remD % p == 0) {
                    remD /= p;
                    count++;
                }
                factors.add(new Factor(p, count));
            }
        }
        if (remD > 1)
            factors.add(new Factor(remD, 1));

        factors.sort((a, b) -> Integer.compare(a.p, b.p));
        List<Factor> merged = new ArrayList<>();
        for (Factor f : factors) {
            if (!merged.isEmpty() && merged.get(merged.size() - 1).p == f.p) {
                merged.get(merged.size() - 1).e += f.e;
            } else {
                merged.add(new Factor(f.p, f.e));
            }
        }
        return merged;
    }

    static void buildDivisorsRecursive(List<Factor> factors, int index, long current, List<Long> out) {
        if (index == factors.size()) {
            out.add(current);
            return;
        }
        long val = current;
        int p = factors.get(index).p;
        int e = factors.get(index).e;
        for (int i = 0; i <= e; ++i) {
            buildDivisorsRecursive(factors, index + 1, val, out);
            val *= p;
        }
    }

    static byte[] squareRootsWithPositiveRep(long limit, int d, int[] spf) {
        int maxRoot = (int) isqrt(limit);
        int maxY = (int) isqrt(limit / d);
        byte[] mark = new byte[maxRoot + 1];

        for (int y = 1; y <= maxY; ++y) {
            long D = (long) d * y * y;
            List<Factor> factors = factorizeDy2(y, d, spf);
            List<Long> divs = new ArrayList<>();
            buildDivisorsRecursive(factors, 0, 1L, divs);

            for (long u : divs) {
                long v = D / u;
                if (u > v)
                    continue;
                if (((u + v) & 1) != 0)
                    continue;
                if (v <= u)
                    continue;

                long m = (u + v) / 2;
                if (m <= maxRoot) {
                    mark[(int) m] = 1;
                }
            }
        }
        return mark;
    }

    static long countGoodSquares(long limit) {
        int maxRoot = (int) isqrt(limit);
        int[] spf = buildSpfTable(maxRoot);

        byte[] d1 = squareRootsWithPositiveRep(limit, 1, spf);
        byte[] d2 = squareRootsWithPositiveRep(limit, 2, spf);
        byte[] d3 = squareRootsWithPositiveRep(limit, 3, spf);
        byte[] d7 = squareRootsWithPositiveRep(limit, 7, spf);

        long count = 0;
        for (int m = 1; m <= maxRoot; ++m) {
            if (d1[m] == 1 && d2[m] == 1 && d3[m] == 1 && d7[m] == 1) {
                count++;
            }
        }
        return count;
    }

    public static String solve() {
        long limit = 2000000000L;
        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());

        long sqKernel = isqrt(limit);
        long singleCont = sumSingletonKernels(limit, threads);
        long multiCont = sumMultiPrimeKernels(limit);

        long intRep = sqKernel + singleCont + multiCont;
        long totalSq = isqrt(limit);
        long goodSq = countGoodSquares(limit);

        return String.valueOf(intRep - totalSq + goodSq);
    }

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