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

Call a positive integer $n$ special if for each of the four coefficients $d=1,2,3,7$ there are positive integers $a_d$ and $b_d$ with $n=a_d^2+d\,b_d^2$; the four representations are independent of one another. The statement gives two examples,

$$\begin{aligned}3600&=48^2+36^2\\&=20^2+2\cdot40^2\\&=30^2+3\cdot30^2\\&=45^2+7\cdot15^2,\\88201&=99^2+280^2\\&=287^2+2\cdot54^2\\&=283^2+3\cdot52^2\\&=197^2+7\cdot84^2,\end{aligned}$$

and says that $75373$ special numbers do not exceed $10^7$. We must count the special numbers that do not exceed $L=2\times10^9$. ^[1]

The whole solution in one sentence: a number that is not a perfect square is special exactly when every prime occurring in it to an odd power leaves the remainder $1$, $25$ or $121$ on division by $168$, so such numbers are counted as $q\,m^2$ with $q$ a product of these primes, while each perfect square $m^2$ is tested separately by factoring $m^2-a^2=(m-a)(m+a)$.

Mathematical Approach

We first decide which numbers have all four representations when zero entries are allowed, and then turn to positive entries, which make a difference only for perfect squares. Throughout, $d$ is one of the coefficients $1,2,3,7$.

1. Small cases

A direct search over all pairs $(a,b)$ finds $34$ special numbers up to $4000$: $26$ primes, from $193$, $337$, $457$, $673$, $1009$ up to $3889$; seven such primes times a square, $772=2^2\cdot193$, $1737=3^2\cdot193$, $3088=4^2\cdot193$, $1348=2^2\cdot337$, $3033=3^2\cdot337$, $1828=2^2\cdot457$ and $2692=2^2\cdot673$; and the square $3600=60^2$. The smallest one is the prime

$$\begin{aligned}193&=12^2+7^2=11^2+2\cdot6^2\\&=1^2+3\cdot8^2=9^2+7\cdot4^2.\end{aligned}$$

The primes leave the remainders $25$, $1$, $121$, $1$, $1$, … on division by $168=2^3\cdot3\cdot7$, and the second example of the statement is $88201=193\cdot457$.

2. When a prime is $a^2+d\,b^2$

Let $p$ be an odd prime that does not divide $d$. Then $p=a^2+d\,b^2$ for some integers $a,b$ if and only if $-d$ is a square modulo $p$, that is, $r^2\equiv-d\pmod p$ for some integer $r$.

If $p=a^2+d\,b^2$, then $p$ does not divide $b$, since otherwise it would divide $a$ too and $p^2$ would divide $p$; so $b$ has an inverse $b^{-1}$ modulo $p$, a number with $b\,b^{-1}\equiv1\pmod p$, and $r=a\,b^{-1}$ works. For the converse we use Thue's lemma ^[2], a pigeonhole argument: two of the $(\lfloor\sqrt p\rfloor+1)^2\gt p$ pairs $(s,t)$ with $0\le s,t\le\lfloor\sqrt p\rfloor$ give the same remainder of $s-rt$ modulo $p$, and their difference is a pair $(x,y)\ne(0,0)$ with $|x|,|y|\lt\sqrt p$ and $x\equiv ry\pmod p$. Then $p$ divides $x^2+d\,y^2\equiv(r^2+d)\,y^2$, and $0\lt x^2+d\,y^2\lt(1+d)\,p$, so $x^2+d\,y^2=kp$ with $1\le k\le d$. For $k=1$ we are done; every other $k$ is impossible or leads back to $p$:

$d=2$, $k=2$: $x$ is even, and $p=y^2+2(x/2)^2$.
$k=d$ for $d=3$ or $7$: $d$ divides $x$, and $p=y^2+d\,(x/d)^2$.
$k$ even for $d=3$ or $7$: $x$ and $y$ have the same parity. For $d=3$ the only even $k$ is $2$; odd $x,y$ give $x^2+3y^2\equiv4\pmod8$, while $2p\equiv2$ or $6$, and even $x,y$ would make $4$ divide $2p$. For $d=7$ and $k=2,4,6$, odd $x,y$ make $x^2+7y^2$ divisible by $8$, which no such $kp$ is; so $x,y$ are even, $4$ divides $kp$, hence $k=4$ and $p=(x/2)^2+7(y/2)^2$.
$d=7$, $k=3$ or $5$: as $-1$ is not a square modulo $3$ and $-2$ is not a square modulo $5$, the congruence $x^2+7y^2\equiv0\pmod k$ forces $k$ to divide $x$ and $y$, so $k^2$ divides $kp$ and $p=k$. But $-7$ is a square modulo neither $3$ nor $5$, so this case does not occur.

For example, $p=11$ and $d=2$ allow $r=3$, as $3^2=9\equiv-2\pmod{11}$. The pairs $(2,0)$ and $(0,3)$ give the same remainder $2$ of $s-3t$; their difference $(x,y)=(2,-3)$ gives $x^2+2y^2=22=2\cdot11$, so $k=2$, and the first case yields $11=3^2+2\cdot1^2$.

3. Which primes qualify

Quadratic reciprocity and its supplementary laws tell when $-d$ is a square modulo an odd prime $p$ not dividing $d$ ^[3]:

$-1$ exactly when $p\equiv1\pmod4$;
$-2$ exactly when $p\equiv1$ or $3\pmod8$;
$-3$ exactly when $p\equiv1\pmod3$;
$-7$ exactly when $p\equiv1$, $2$ or $4\pmod7$, the nonzero squares modulo $7$.

With Section 2 these are the classical descriptions of the primes $a^2+b^2$, $a^2+2b^2$ and $a^2+3b^2$, announced by Fermat ^[4], and of the primes $a^2+7b^2$ ^[5]; for example, $11\equiv4\pmod7$ and $11=2^2+7\cdot1^2$.

4. Good primes

A prime has all four representations exactly when it leaves the remainder $1$, $25$ or $121$ on division by $168$; we call such primes good. The primes $2$, $3$, $7$ fail: $2$ is not $a^2+3b^2$, and $3$ and $7$ are not sums of two squares. Any other prime $p$ divides none of the coefficients, so by Sections 2 and 3 it qualifies exactly when $p\equiv1\pmod8$, the only class allowed for both $-1$ and $-2$, when $p\equiv1\pmod3$ and when $p\equiv1$, $2$ or $4\pmod7$. By the Chinese remainder theorem ^[6] the first two give $p\equiv1\pmod{24}$, and together with the third exactly the remainders $1$, $25=1+24$ and $121=1+5\cdot24$ modulo $168$, which are $1$, $4$ and $2$ modulo $7$. They are the squares $1^2$, $5^2$, $11^2$, the only squares modulo $168$ of numbers prime to $168$, and the good primes are also the primes $x^2+168y^2$ ^[7]. The numbers $1$, $25$, $121$ and $169=13^2$ are not prime, so $193$ is the smallest good prime.

5. Multiplying representations

Brahmagupta's identity ^[8]

$$\begin{aligned}&(a^2+d\,b^2)(c^2+d\,e^2)\\&\quad=(ac-d\,be)^2+d\,(ae+bc)^2\\&\quad=(ac+d\,be)^2+d\,(ae-bc)^2\end{aligned}$$

shows that a product of numbers of the form $a^2+d\,b^2$ has the same form, and so has such a number times a square: $n\,m^2=(am)^2+d\,(bm)^2$. With

$$\begin{aligned}457&=21^2+4^2=13^2+2\cdot12^2\\&=5^2+3\cdot12^2=3^2+7\cdot8^2,\end{aligned}$$

the identity, with $(a,b)$ taken from $193$ and $(c,e)$ from $457$, turns their representations into exactly the four representations of $88201$ in the statement. The last line gives $12\cdot21+7\cdot4=280$ and $12\cdot4-7\cdot21=-99$ for $d=1$, and $11\cdot13+2\cdot6\cdot12=287$ and $11\cdot12-6\cdot13=54$ for $d=2$; the middle line gives $1\cdot5-3\cdot8\cdot12=-283$ and $1\cdot12+8\cdot5=52$ for $d=3$, and $9\cdot3-7\cdot4\cdot8=-197$ and $9\cdot8+4\cdot3=84$ for $d=7$.

6. Primes to an odd power

Now let $n=a^2+d\,b^2$ with integers $a,b\ge0$, and let $p$ be an odd prime not dividing $d$ for which $-d$ is not a square modulo $p$. If $p$ divides $n$, then $p$ divides $b$, as otherwise $r=a\,b^{-1}$ would satisfy $r^2\equiv-d$; hence $p$ divides $a$, and $n/p^2=(a/p)^2+d\,(b/p)^2$ has the same form. Each repetition of this descent removes $p^2$, so $p$ occurs in $n$ to an even power. The prime $2$ also occurs in $a^2+3b^2$ to an even power: if $a$ and $b$ are even we divide by $4$, if exactly one is odd then $n$ is odd, and if both are odd then $a^2+3b^2\equiv4\pmod8$, so $2$ occurs exactly twice.

So if $n$ has all four representations, every prime $p$ occurring in $n$ to an odd power is good: $p\ne2$ by the case $d=3$; $p\ne3,7$ by the case $d=1$, as $-1$ is not a square modulo $3$ or $7$; and any other $p$ divides no coefficient, so $-d$ is a square modulo $p$ for all four $d$, and $p$ is good by Sections 3 and 4.

7. The squarefree part decides

Every positive integer is, in exactly one way, $n=q\,m^2$ with $q$ squarefree, that is, not divisible by the square of a prime ^[9]; $q$ is the product of the primes that occur in $n$ to an odd power, for example $3088=193\cdot4^2$ and $3600=1\cdot60^2$. The code calls $q$ the kernel of $n$.

Theorem. The number $n=q\,m^2$ has, for each $d$, a representation $n=a^2+d\,b^2$ with integers $a,b\ge0$ if and only if every prime factor of $q$ is good.

Section 6 proves "only if". Conversely, every good prime has all four representations by Section 4, so their product $q$ has them by Section 5 (for $q=1$ take $1=1^2+d\cdot0^2$), and so does $q\,m^2$.

8. Positive entries

Let $n=q\,m^2$ satisfy the theorem. If $q\gt1$, then $n$ is not a square, so no representation has $b=0$; and $a=0$ would mean $n=d\,b^2$, a square for $d=1$, while for $d=2,3,7$ the prime $d$ would divide $q$, although $2$, $3$, $7$ are not good. Hence every such $n$ with $q\gt1$ is special. A square, however, always has the four representations $m^2=m^2+d\cdot0^2$ with a zero entry, so squares need their own test. With $A(L)$ the number of special numbers up to $L$ and $G(L)$ the number of special squares among them,

$$A(L)=\sum_{q}\left\lfloor\sqrt{L/q}\right\rfloor+G(L),$$

where $q\gt1$ runs over the squarefree numbers made of good primes: for fixed $q$ the numbers $q\,m^2\le L$ are those with $m\le\sqrt{L/q}$, and by uniqueness different $q$ give different numbers.

9. Kernels with one prime and with several

Most terms have $q=p$, a single good prime $p\le L$. In a product of two or more good primes every factor is at most $L/193$, since the others contribute at least $193$. The smallest such products are $65041=193\cdot337$ and $88201=193\cdot457$. Up to $L$ there are products of two and of three good primes but none of four, as $193\cdot337\cdot457=29723737$ is below $L$ and $193\cdot337\cdot457\cdot673=20004075001$ is above it.

10. Squares and factor pairs

For a square, $m^2=a^2+d\,b^2$ with $a,b\gt0$ means $(m-a)(m+a)=d\,b^2$. With $u=m-a$ and $v=m+a$ this is a factorization $d\,b^2=uv$ with $0\lt u\lt v$, and $u\equiv v\pmod2$ because $u+v=2m$. Conversely, each such factorization gives positive integers $m=(u+v)/2$ and $a=(v-u)/2$ with $m^2-a^2=d\,b^2$. So we find all $m$ with a positive representation of $m^2$ for a given $d$ by running through $b\le\sqrt{L/d}$, listing the divisors $u$ of $d\,b^2$ and keeping the pairs $u\lt v$ of equal parity. For $m=60$:

$d=1$, $b=36$: $1296=12\cdot108$, so $a=48$;
$d=2$, $b=40$: $3200=40\cdot80$, so $a=20$;
$d=3$, $b=30$: $2700=30\cdot90$, so $a=30$;
$d=7$, $b=15$: $1575=15\cdot105$, so $a=45$.

These are the representations of $3600$ in the statement. The smallest special squares are $60^2$, $68^2$, $110^2$, $111^2$, $120^2$, and $139$ of the $1000$ squares up to $10^6$ are special.

11. Worked example: $L=4000$

There are $26$ good primes up to $4000$. The terms $\lfloor\sqrt{4000/p}\rfloor$ are $4$ for $p=193$, counting $193$, $772$, $1737$, $3088$; $3$ for $p=337$; $2$ for $p=457$ and $p=673$; and $1$ for the other $22$, as $4\cdot1009=4036\gt4000$. No product of two good primes is that small, so the sum is $4+3+2+2+22=33$. Only $60^2$ is a special square, so $G(4000)=1$ and $A(4000)=34$, as found in Section 1.

12. The checkpoint

For $L=10^7$: $41310$ good primes contribute $73297$; $1106$ products of two good primes contribute $1467$, while three good primes already have a product above $10^7$; and $609$ of the $3162$ squares are special. Together $73297+1467+609=75373$, the value in the statement.

How the Code Works

The C++ program ^[10] starts in main, which reads the optional arguments --limit=N and --skip-checkpoints, runs run_checkpoints unless told to skip them, and prints the result of count_special_numbers. That function computes Section 8 as $R-\lfloor\sqrt L\rfloor+G(L)$: count_integer_representable counts the $R$ numbers $n\le L$ that satisfy the theorem of Section 7, the squares included through the term $\lfloor\sqrt L\rfloor$ for $q=1$; the squares are subtracted again, and count_good_squares adds back the special ones. isqrt_u64 corrects a floating-point square root up or down; taking it of $\lfloor L/q\rfloor$ instead of $L/q$ is safe, since $k^2\le\lfloor x\rfloor$ and $k^2\le x$ agree for integers $k$.

sum_singleton_kernels finds the base primes up to $\lfloor\sqrt L\rfloor$ with sieve_primes_up_to, cuts the range from $2$ to $L$ into equal blocks, and lets singleton_sum_in_range sieve each block in segments of $2^{20}$ numbers ^[11], each base prime $p$ crossing out its multiples from $p^2$ or from its first multiple in the segment, whichever is larger. For every surviving prime, is_good_prime checks $p\gt7$ and the remainder modulo $168$, and a good prime adds $\lfloor\sqrt{L/p}\rfloor$.

sum_multi_prime_kernels returns $0$ if $L\lt193^2$; otherwise collect_good_primes_up_to sieves up to $\lfloor L/193\rfloor$ and keeps the good primes, and sum_multi_prime_kernels_dfs builds products of distinct good primes by depth-first search, always appending a larger prime, so that each product is reached once. Each loop stops at the first $p$ for which the current product exceeds $\lfloor L/p\rfloor$, a test that cannot overflow, and each product $q$ of at least two primes adds $\lfloor\sqrt{L/q}\rfloor$.

count_good_squares implements Section 10. build_spf_table stores the smallest prime factor of every number up to $\lfloor\sqrt L\rfloor$. For each $d$, square_roots_with_positive_rep runs through $b\le\lfloor\sqrt{L/d}\rfloor$, named y in the code: factorize_dy2 factors $d\,b^2$ by doubling the exponents of $b$ and adding the primes of $d$, build_divisors_recursive lists the divisors $u$, and each pair $u\lt v=d\,b^2/u$ of equal parity marks $m=(u+v)/2$ if $m\le\lfloor\sqrt L\rfloor$. The $m$ marked for all four $d$ are counted.

run_checkpoints checks with has_positive_rep_bruteforce that $3600$ and $88201$ have all four positive representations, compares the fast count for $L=200000$ with brute_force_count, which marks every $a^2+d\,b^2\le L$ with $a,b\ge1$, and checks $75373$ for $L=10^7$; a failure stops the program. The Python and Java versions mirror these functions, Java in camel case and Python with the search nested in solve, but compute only $L=2\times10^9$ and run no checkpoints. All three use exact integers: unsigned 64-bit in C++ and long in Java for products and sums, unbounded integers with math.isqrt in Python. Kernels, values $d\,b^2$ and their divisors never exceed $L$; the largest intermediate numbers arise in the recursive divisor builders of C++ and Java, which multiply once more after the last power of each prime, and they stay below $L\sqrt L\lt9\times10^{13}$, far inside 64 bits.

Complexity and Verification

The one-prime part dominates. Its segmented sieve crosses out the multiples of the $4648$ primes up to $\lfloor\sqrt L\rfloor=44721$, which takes $O(L\log\log L)$ operations ^[11], and scans every number from $2$ to $L$; since a segment of $2^{20}$ numbers is longer than $\sqrt L$, going through the base primes once per segment adds only $O(L)$ work. Of its $98222287$ primes, $6136501$ are good, about one in sixteen, just as $3$ of the $48$ remainders prime to $168$ are good; each good prime needs one square root. The second sieve runs up to $\lfloor L/193\rfloor=10362694$ with about $10$ MB of flags, finds $687054$ primes and keeps $42705$ good ones; the search then visits each of the $330749$ products of two and $1337$ products of three good primes once. The square part runs through $44721+31622+25819$${}+16903=119065$ values of $b$ and generates $7033639$ divisors, at most $1080$ for one $d\,b^2$. In total the time is $O(L\log\log L)$. The largest data belong to the second sieve: its flags and its list of primes, which in Python consists of separate integer objects and takes about $25$ MB. In a test on a desktop computer the optimized C++ program needed less than four seconds, checkpoints included.

Separate programs written for this article checked the descriptions of Section 3 for all primes below $10^6$; the cases of Section 2, by listing all pairs $(x,y)$ of Thue's lemma for the primes below $3000$, which give only $k=1,2$ for $d=2$, $k=1,3$ for $d=3$ and $k=1,4,7$ for $d=7$; the theorem of Section 7 and the positivity argument for every $n\le10^6$; and, with special squares found by direct search instead of factor pairs, $A(L)=15$, $34$, $196$, $956$, $1878$, $8637$, $75373$ for $L=2000$, $4000$, $20000$, $10^5$, $2\cdot10^5$, $10^6$, $10^7$, the same as a brute force over all pairs $(a,b)$. At the target, each of the three parts of the C++ code agrees with an independent computation, and the three programs print the same result.

Footnotes and References

The references below supply the classical facts used above. The case analysis after Thue's lemma, the reduction to squarefree parts, the positivity argument and the factor-pair test for squares are derived in this article.

Project Euler 229 — Four Representations Using Squares. The official statement asks how many n ≤ 2 × 10⁹ are a₁² + b₁², a₂² + 2b₂², a₃² + 3b₃² and a₇² + 7b₇² with positive aₖ, bₖ, with the examples 3600 and 88201 and the count 75373 up to 10⁷.
Wikipedia — Thue's lemma. For integers a and m > 1 and positive X, Y with X ≤ m < XY, there are integers x, y with ay ≡ x (mod m), |x| < X and 0 < y < Y. Thue's proof of 1902 is a pigeonhole argument.
Wikipedia — Quadratic reciprocity. For an odd prime p, −1 is a quadratic residue exactly when p ≡ 1 (mod 4), −2 when p ≡ 1 or 3 (mod 8), −3 when p ≡ 1 (mod 3), and −7 when p is a residue modulo 7.
Wikipedia — Fermat's theorem on sums of two squares. An odd prime is a sum of two squares exactly when it is ≡ 1 (mod 4); for odd primes p, Fermat also stated in a letter to Pascal of 1654 that p = x² + 2y² exactly when p ≡ 1 or 3 (mod 8), and p = x² + 3y² exactly when p ≡ 1 (mod 3), where p = 3 = 0² + 3·1² must be set aside.
OEIS A033207 — Primes of the form x² + 7y². They are the primes congruent to 1, 7, 9, 11, 15, 23 or 25 modulo 28, that is 7 and the odd primes ≡ 1, 2 or 4 (mod 7); the entry refers to D. A. Cox, Primes of the Form x² + ny² (Wiley, 1989).
Wikipedia — Chinese remainder theorem. The remainders of an integer modulo pairwise coprime moduli determine its remainder modulo their product uniquely.
OEIS A139506 — Primes of the form x² + 26xy + y². The entry, 193, 337, 457, 673, 1009, …, notes that these are also the primes x² + 168y², and that they are the primes congruent to 1, 25 or 121 modulo 168.
Wikipedia — Brahmagupta's identity. (a² + nb²)(c² + nd²) = (ac ± nbd)² + n(ad ∓ bc)², with matching signs, so products of numbers a² + nb² keep that form; for n = 1 it is the Brahmagupta–Fibonacci identity.
Wikipedia — Square-free integer. Each positive integer is in a unique way the product of a largest square and a square-free integer, its square-free part.
C++ — Euler229.cpp. The linked immutable revision contains the C++ functions described above, from count_special_numbers to run_checkpoints; it is the implementation record for this article.
Wikipedia — Sieve of Eratosthenes. Crossing out the multiples of p may start at p², the sieve takes O(n log log n) operations, and the segmented sieve treats the range one segment at a time with the primes up to the square root.

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());
    }
}