Problem 233: Lattice Points on a Circle

Project Euler Problem 233 Solution

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

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

Detailed mathematical approach

Problem Summary

For a positive integer $N$, the points $(0,0)$, $(N,0)$, $(0,N)$ and $(N,N)$ are the corners of a square, and exactly one circle passes through all four of them. Let $f(N)$ be the number of points with integer coordinates, called lattice points, on this circle. The statement gives $f(10000)=36$ and asks for the sum of all positive integers $N\le10^{11}$ with $f(N)=420$. ^[1] We write $S(L)$ for the sum of all positive integers $N\le L$ with $f(N)=420$, so the task is to compute $S(10^{11})$.

The whole solution in one sentence: the lattice points correspond to the ways of writing $N^2$ as a sum of two squares, their number depends only on the exponents of the primes $p\equiv1\pmod4$ in $N$, and so we list every product of such prime powers up to $10^{11}$ whose exponents give $420$ and multiply it by the sum of its admissible cofactors, read from a single table of running sums.

Mathematical Approach

We turn the lattice points into sums of two squares, count them with Gaussian integers, find every exponent pattern that gives $420$, and add up the admissible numbers grouped by the part of $N$ that decides $f(N)$.

1. The circle and a small example

The centre of the circle is equally far from all four corners, so it is the centre $(N/2,N/2)$ of the square, and the radius is half the diagonal, $N/\sqrt2$. The circle is therefore

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

Multiplying out, the two terms $N^2/4$ on the left cancel against the right-hand side, and the equation becomes

$$x^2+y^2=N(x+y).$$

For $N=5$ a short search finds $12$ lattice points: the four corners and the eight points $(6,2)$, $(6,3)$, $(2,6)$, $(3,6)$, $(-1,2)$, $(-1,3)$, $(2,-1)$ and $(3,-1)$. For example, $(6,2)$ gives $36+4=40=5\cdot8$, and $(-1,2)$ gives $1+4=5=5\cdot1$. So $f(5)=12$.

2. From the circle to sums of two squares

Put $a=x+y-N$ and $b=x-y$. Expanding the squares gives

$$a^2+b^2-N^2=2\bigl(x^2+y^2-N(x+y)\bigr),$$

so a point lies on the circle exactly when $a^2+b^2=N^2$. Conversely, integers $a,b$ with $a^2+b^2=N^2$ give back $x=(N+a+b)/2$ and $y=(N+a-b)/2$, and these are integers: every integer $k$ satisfies $k^2\equiv k\pmod2$, that is, $k^2$ and $k$ leave the same remainder on division by $2$, so $a+b\equiv a^2+b^2=N^2\equiv N\pmod2$. Hence $N+a+b$ is even, and so is $N+a-b$, which differs from it by $2b$. The two substitutions undo each other, so

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

where $r_2(n)$ is the number of pairs of integers $(a,b)$ with $a^2+b^2=n$, counting order and signs. ^[3] Geometrically, the substitution is a similarity with ratio $\sqrt2$ that moves the centre to the origin, so $f(N)$ is also the number of lattice points on the circle of radius $N$ around the origin. ^[4] For $N=5$ the twelve pairs are $(\pm5,0)$, $(0,\pm5)$, $(\pm3,\pm4)$ and $(\pm4,\pm3)$; for example, $(a,b)=(3,4)$ corresponds to $(x,y)=(6,2)$.

3. Gaussian integers

To count these pairs we use the Gaussian integers, the numbers $a+bi$ with integers $a,b$ and $i^2=-1$. The norm of $z=a+bi$ is $z\bar z=a^2+b^2$, where $\bar z=a-bi$ is the conjugate, and the norm of a product is the product of the norms. The pairs counted by $r_2(N^2)$ are therefore exactly the Gaussian integers of norm $N^2$. A unit is a Gaussian integer whose inverse is a Gaussian integer too, and a Gaussian prime is a nonzero non-unit that is not a product of two non-units. In the same way, $p\equiv1\pmod4$ and $q\equiv3\pmod4$ mean that $p$ and $q$ leave the remainders $1$ and $3$ on division by $4$. Three facts do all the work: ^[2]

the units, the Gaussian integers of norm $1$, are $1$, $-1$, $i$ and $-i$;
every nonzero Gaussian integer is a product of Gaussian primes in only one way, apart from the order of the factors and from unit factors;
apart from unit factors, the Gaussian primes are $1+i$, of norm $2$; every ordinary prime $q\equiv3\pmod4$, of norm $q^2$; and, for every ordinary prime $p\equiv1\pmod4$, two primes $\pi$ and $\bar\pi$ with $\pi\bar\pi=p$ that are not unit multiples of each other.

For example, $5=(2+i)(2-i)$ and $13=(3+2i)(3-2i)$ split, $3$ does not, and $2=-i(1+i)^2$.

4. Counting the lattice points

Write the prime factorisation of $N$ as

$$N=2^{\alpha}\,p_1^{e_1}\cdots p_k^{e_k}\,q_1^{g_1}\cdots q_l^{g_l}$$

with every $p_j\equiv1\pmod4$ and every $q_j\equiv3\pmod4$, and split $p_j=\pi_j\bar\pi_j$ as in Section 3. Claim:

$$f(N)=4\,(2e_1+1)(2e_2+1)\cdots(2e_k+1).$$

Proof. Let $z$ be a Gaussian integer of norm $N^2$. Since $z$ divides $z\bar z=N^2$, its prime factors are among those of $N^2$, namely $1+i$, the $\pi_j$, the $\bar\pi_j$ and the $q_j$. So

$$z=u\,(1+i)^{s}\prod_j\pi_j^{x_j}\,\bar\pi_j^{\,y_j}\prod_j q_j^{t_j}$$

with a unit $u$ and exponents $s,x_j,y_j,t_j\ge0$. Taking norms gives

$$2^{s}\prod_j p_j^{x_j+y_j}\prod_j q_j^{2t_j}=N^2,$$

and comparing the exponent of each ordinary prime with its exponent in $N^2$ gives $s=2\alpha$, $t_j=g_j$ and $x_j+y_j=2e_j$. Conversely, with $s=2\alpha$ and $t_j=g_j$, every unit $u$ and every choice of $x_j\in\{0,1,\dots,2e_j\}$ with $y_j=2e_j-x_j$ gives a Gaussian integer of norm $N^2$, and different choices give different numbers, again by unique factorisation, because $\pi_j$ and $\bar\pi_j$ are not unit multiples of each other. There are $4$ units and $2e_j+1$ choices for each $p_j$, while the factors $2$ and $q_j$ leave no choice, which proves the claim. If $N$ has no prime $p_j$, the product is empty and equals $1$, so $f(N)=4$.

For $N=5$ the twelve numbers are $u(2+i)^2=u(3+4i)$, $u(2+i)(2-i)=5u$ and $u(2-i)^2=u(3-4i)$ with the four units $u$, which are the twelve pairs of Section 2. The checkpoint of the statement agrees: $10000=2^4\cdot5^4$, so $f(10000)=4\cdot(2\cdot4+1)=36$. A prime $q\equiv3\pmod4$ adds nothing: $f(3)=4$, from $(\pm3,0)$ and $(0,\pm3)$ alone. The claim is the product formula for $r_2(n)$ with $n=N^2$. ^[3]

5. Right triangles with hypotenuse $N$

Apart from $(\pm N,0)$ and $(0,\pm N)$, every solution of $a^2+b^2=N^2$ has $a,b\ne0$ and $|a|\ne|b|$, since $2a^2=N^2$ has no solution in positive integers ($\sqrt2$ is irrational). Such solutions come in groups of eight, $(\pm a,\pm b)$ and $(\pm b,\pm a)$, one group for every right triangle with integer legs $0\lt a\lt b$ and hypotenuse $N$. Hence $f(N)=4+8H(N)$, where $H(N)$ is the number of these triangles, and Section 4 gives $H(N)=\bigl(\prod_j(2e_j+1)-1\bigr)/2$, the formula of OEIS A046080. ^[5] For $N=10000$ there are $H=4$ triangles, with legs $(2800,9600)$, $(3520,9360)$, $(5376,8432)$ and $(6000,8000)$. Since $420=4+8\cdot52$, the problem asks for the hypotenuses $N\le10^{11}$ of exactly $52$ right triangles with integer sides.

6. Exponent patterns for 420

By Section 4, $f(N)=420$ means

$$(2e_1+1)(2e_2+1)\cdots(2e_k+1)=105=3\cdot5\cdot7.$$

Each factor is odd and at least $3$, and conversely every divisor greater than $1$ of the odd number $105$ has this form, so we need every way of writing $105$ as a product of factors greater than $1$, ignoring their order. Each factor collects some of the three primes $3$, $5$ and $7$, so the ways correspond to the groupings of these primes: one factor, $105$; two factors, $3\cdot35$, $5\cdot21$ or $7\cdot15$; or three factors, $3\cdot5\cdot7$. Solving $2e+1=d$ for every factor $d$ gives the five exponent patterns

$$(52),\quad(17,1),\quad(10,2),\quad(7,3),\quad(3,2,1).$$

So $f(N)=420$ holds exactly when the primes $p\equiv1\pmod4$ of $N$ carry one of these lists of exponents, in some order, while the exponent of $2$ and those of the primes $q\equiv3\pmod4$ are arbitrary.

7. Which patterns fit below $10^{11}$

The smallest number with a given pattern uses the smallest primes $\equiv1\pmod4$, which are $5$, $13$ and $17$, and puts the largest exponent on the smallest of them. Indeed, if $p\lt q$ and $e\lt e'$, then $p^{e'}q^{e}\lt p^{e}q^{e'}$, because the quotient of the two sides is $(p/q)^{e'-e}\lt1$. The smallest values are

$(52)$: $5^{52}\approx2.2\cdot10^{36}$;
$(17,1)$: $5^{17}\cdot13=9\,918\,212\,890\,625$;
$(10,2)$: $5^{10}\cdot13^2=1\,650\,390\,625$;
$(7,3)$: $5^{7}\cdot13^3=171\,640\,625$;
$(3,2,1)$: $5^{3}\cdot13^2\cdot17=359\,125$.

The first two exceed $10^{11}$, so only the last three patterns occur, and $359\,125$ is the smallest $N$ with $f(N)=420$.

8. Cores and allowed multipliers

Call a positive integer an allowed multiplier if it has no prime factor $\equiv1\pmod4$; the first ones are $1,2,3,4,6,7,8,9,11,12,14,16$. By Sections 4 and 5 these are exactly the numbers $m$ with $f(m)=4$, which are the hypotenuse of no right triangle with integer sides, the nonhypotenuse numbers of OEIS A004144. ^[6] Every positive integer $N$ splits in exactly one way as $N=c\,m$, where the core $c$ is the product of the full prime powers in $N$ whose primes are $\equiv1\pmod4$, and $m$ is an allowed multiplier. The formula of Section 4 only sees the core, so $f(N)=f(c)$.

Let $\mathcal C$ be the set of cores whose exponents form one of the patterns of Section 6, and let $A(X)$ be the sum of all allowed multipliers $m\le X$. The numbers $N\le L$ with $f(N)=420$ are exactly the products $c\,m$ with $c\in\mathcal C$ and an allowed multiplier $m\le\lfloor L/c\rfloor$, and each of them arises from only one such pair. Adding them core by core gives

$$S(L)=\sum_{c\in\mathcal C,\ c\le L}c\cdot A\bigl(\lfloor L/c\rfloor\bigr).$$

The first values are $A(1)=1$, $A(2)=3$, $A(3)=6$, $A(4)=A(5)=10$ and $A(6)=16$; $5$ is not allowed.

9. Worked example: $L=10^6$

Six cores lie below $10^6$:

$359125=5^3\cdot13^2\cdot17$ and $469625=5^3\cdot13\cdot17^2$, with $\lfloor10^6/c\rfloor=2$;
$612625=5^3\cdot13^2\cdot29$, $781625=5^3\cdot13^2\cdot37$, $866125=5^3\cdot13^2\cdot41$ and $933725=5^2\cdot13^3\cdot17$, with $\lfloor10^6/c\rfloor=1$.

With $A(2)=3$ and $A(1)=1$, the formula of Section 8 gives

$$\begin{aligned}S(10^6)&=3\,(359125+469625)\\&\quad+612625+781625+866125+933725\\&=5\,680\,350.\end{aligned}$$

The eight numbers behind this sum are the six cores, $718250=2\cdot359125$ and $939250=2\cdot469625$; the next candidate, $3\cdot359125$, already exceeds $10^6$. Computing $f(N)$ directly for every $N\le10^6$ finds exactly these eight numbers.

10. How far the search must go

In a core of pattern $(3,2,1)$, the two primes with exponents $3$ and $2$ contribute at least $5^3\cdot13^2=21125$: the other order gives $13^3\cdot5^2=54925$, and without the prime $5$ both primes are at least $13$. So the prime with exponent $1$ is at most $L/21125$, which is about $4.73\cdot10^6$ for $L=10^{11}$. Every prime with exponent at least $2$, in any pattern, satisfies $p^2\le L$, so for $L=10^{11}$ it is at most $316\,227$, far below this bound. By Section 7 every core is at least $359125$, so no multiplier larger than $M=\lfloor L/359125\rfloor$ is ever needed; for $L=10^{11}$, $M=278\,454$.

For $L=10^{11}$ there are $585\,259$ cores: $585\,235$ of pattern $(3,2,1)$, $13$ of pattern $(7,3)$ and $11$ of pattern $(10,2)$. The largest prime that occurs is $4\,733\,689$, in the core $5^3\cdot13^2\cdot4733689$.

11. Listing the cores

A core is built by choosing primes $p_1\lt p_2\lt\cdots$ among the primes $\equiv1\pmod4$ and giving them the exponents of a pattern in a fixed order. Running through all orders of all patterns produces every core exactly once, because the factorisation of a core fixes its increasing list of primes and the order of their exponents. The primes are chosen one at a time, from small to large, and two tests end the loop over the current prime $p$ early:

$p^e$ exceeds the room $\lfloor L/P\rfloor$ left by the product $P$ of the earlier prime powers;
the smallest completion, which gives the remaining exponents to the next primes in increasing order, does not fit into the room left after $p^e$.

A larger $p$ makes $p^e$, the product and the smallest completion larger, so a test that fails keeps failing for all later primes, and stopping the loop loses no core. Every branch that passes both tests ends in at least one core, so the work is proportional to the number of cores: for $L=10^{11}$ the search function is called $585\,830$ times, and $585\,259$ of these calls end in a core.

12. Running sums of the allowed multipliers

To decide which $m\le M$ are allowed, the program uses the smallest prime factor $\operatorname{spf}(n)$ of every $n\le M$. The linear sieve computes all of them in $O(M)$ steps: it runs through $n=2,3,\dots,M$, records $n$ as a prime if nothing has marked it yet, and then sets $\operatorname{spf}(p\,n)=p$ for every recorded prime $p\le\operatorname{spf}(n)$ with $p\,n\le M$. Every composite number is marked exactly once, as its smallest prime factor times a number with no smaller prime factor. ^[8] The prime factors of $n\gt1$ are $\operatorname{spf}(n)$ and those of $n/\operatorname{spf}(n)$, so $n$ is allowed exactly when $\operatorname{spf}(n)$ is $2$ or $\equiv3\pmod4$ and $n/\operatorname{spf}(n)$ is allowed; starting from the allowed number $1$, this decides every $n$ in increasing order. For example, $12$ is allowed because $\operatorname{spf}(12)=2$ and $6$ is allowed, while $15$ is not: $\operatorname{spf}(15)=3$ is fine, but $15/3=5$ is not allowed. A running sum then gives $A(X)$ for every $X\le M$. For $L=10^{11}$ there are $82\,107$ allowed multipliers up to $278\,454$, with sum $A(278454)=11\,144\,226\,412$.

How the Code Works

The C++ program ^[9] accepts --limit= to set $L$, by default $10^{11}$, and --skip-checkpoints. Its function solve_sum carries out Sections 8 to 12:

estimate_prime_bound returns $\max(200,\lfloor L/21125\rfloor+100)$, the bound of Section 10 with a margin, and primes_1mod4_up_to keeps the primes $\equiv1\pmod4$ produced by sieve_primes, a sieve of Eratosthenes that crosses out the multiples of each prime $p$ from $p^2$ on. ^[7] For $L=10^{11}$ these are $165\,553$ of the $331\,181$ primes up to $4\,733\,827$.
enumerate_core_numbers runs through the orders of the five patterns of Section 6 with std::next_permutation and calls enumerate_core_for_order, the search of Section 11. pow_with_limit computes $p^e$ but returns the limit plus one as soon as the power would exceed it, so no product overflows, and tail_fits tests the smallest completion. The patterns $(52)$ and $(17,1)$ fail at their first prime. The cores are then sorted and duplicates are removed, a safeguard, since by Section 11 there are none.
build_allowed_multiplier_prefix receives $M=\lfloor L/c\rfloor$ for the smallest core $c$, builds the table of smallest prime factors with build_spf, the linear sieve of Section 12, fills the byte array valid by the rule of Section 12 and stores $A(0),\dots,A(M)$ as unsigned 128-bit integers, a GCC and Clang extension. ^[10]
accumulate_range adds $c\cdot A(\lfloor L/c\rfloor)$ over a range of cores. By default the list of cores is cut into consecutive blocks whose sums are added at the end; every core lies in exactly one block, so the cut does not change the total.

The main program prints the 128-bit total, converted to decimal by to_string_u128, and the elapsed time. Before that, unless the checkpoints are skipped, run_checkpoints checks four things. f_value_from_factorization, which factors $N$ by trial division and multiplies by $2e+1$ for each prime $\equiv1\pmod4$ that divides $N$ exactly $e$ times, must give $f(10000)=36$, and for every $N\le400$ it must agree with brute_circle_count, which counts the solutions of $a^2+b^2=N^2$ directly; this tests Section 4 and takes Section 2 for granted. brute_sum_direct factors every $n\le2\cdot10^6$ with a table of smallest prime factors and adds the $n$ with $f(n)=420$, and its sum must equal that of solve_sum; both give $S(2\cdot10^6)=37\,661\,475$, from $29$ numbers. Finally solve_sum computes $S(10^{10})$ in one block and with the default blocks, and the two results must agree; both are $2\,709\,839\,225\,627\,731$. With a smaller --limit=, these two checks use $\min(L,2\cdot10^6)$ and $\min(L,10^{10})$ instead.

The Python and Java programs contain the same main functions, which take the limit as a parameter, but their entry points compute only $L=10^{11}$ and run no checkpoints. Python's integers have no fixed size, and it takes the orders of a pattern from itertools.permutations. Java writes its own nextPermutation, keeps primes, smallest prime factors and flags in int and byte arrays, with room in buildSpf for $\lfloor M/10\rfloor+100=27\,945$ primes, more than the $24\,311$ primes up to $M$, and computes the products and sums in long, whose operations do not report overflow. ^[10] This is safe: $A(M)=11\,144\,226\,412$, every term $c\cdot A(\lfloor L/c\rfloor)$ is below $4.1\cdot10^{15}$, and the total is at most $L$ times the number of products $c\,m$, which is at most $\sum_c\lfloor L/c\rfloor=11\,218\,803$. So every partial sum stays below $1.2\cdot10^{18}$, far below $2^{63}-1\approx9.22\cdot10^{18}$.

Complexity and Verification

Let $B\approx L/21125$ be the prime bound, $C$ the number of cores and $M=\lfloor L/359125\rfloor$. The sieve of Eratosthenes takes $O(B\log\log B)$ time and $B$ bytes, ^[7] the search of Section 11 takes $O(C)$ steps and sorting its result $O(C\log C)$, the linear sieve and the running sums take $O(M)$, and the final sum takes one table lookup per core. For $L=10^{11}$ we have $B=4\,733\,827$, $C=585\,259$ and $M=278\,454$: instead of examining $10^{11}$ numbers, the program handles a few million table entries. In the main computation of the C++ program the sieve, the cores and the 128-bit running sums hold about $4.7$, $4.7$ and $4.5$ MB of data; the default checkpoints add a table of smallest prime factors for $2\cdot10^6$ numbers, about $8$ MB. On the Apple M4 Max used for this article the C++ program needs $0.03$ seconds without checkpoints and $0.07$ seconds with them, the Java program about $0.1$ seconds and the Python program under one second.

All three programs give the same result, and the C++ checkpoints pass. For this article we also checked, with programs that are not part of the published solutions:

a direct count of the lattice points on $x^2+y^2=N(x+y)$ agrees with Section 4 for every $N\le400$ and for $N=10000$;
a separate search without the completion test finds the same $585\,259$ cores, and the published search produces none of them twice;
crossing out the multiples of all primes $\equiv1\pmod4$ gives the same allowed multipliers as the rule of Section 12;
computing $f(N)$ for every $N\le L$ agrees with the core method and with the C++ program run with --limit= for $L=10^6$, $2\cdot10^6$, $10^7$, $10^8$ and $2\cdot10^9$, for instance $S(10^7)=1\,722\,531\,830$ and $S(2\cdot10^9)=106\,536\,068\,476\,048$; the last limit already contains cores of all three patterns, since $5^{10}\cdot13^2\lt2\cdot10^9$;
for $L=10^{11}$, the separate search with the separate allowed multipliers gives the same value as the three programs.

Footnotes and References

The references below supply the background facts used above. The substitution of Section 2, the count of Section 4, the exponent patterns, the split into cores and allowed multipliers and the bounds of the search are derived in this article.

Project Euler 233 — Lattice Points on a Circle. The official statement defines f(N) as the number of points with integer coordinates on the circle through (0, 0), (N, 0), (0, N) and (N, N), gives f(10000) = 36 and asks for the sum of all positive integers N ≤ 10¹¹ with f(N) = 420.
Wikipedia — Gaussian integer. The norm a² + b² of a + bi is multiplicative, the units are 1, −1, i and −i, and factorisation into Gaussian primes is unique. 2 = (1 + i)(1 − i) ramifies, a prime p ≡ 3 (mod 4) stays prime, and a prime p ≡ 1 (mod 4) is the product of two conjugate Gaussian primes that are not associates, as in 5 = (2 + i)(2 − i).
Wikipedia — Sum of squares function. r₂(n) counts the representations n = a² + b² with order and signs. The article gives r₂(n) = 4(d₁(n) − d₃(n)) and r₂(n) = 4(h₁ + 1)(h₂ + 1)⋯ over the exponents hᵢ of the primes ≡ 1 (mod 4) in n, valid when every prime ≡ 3 (mod 4) has an even exponent; otherwise r₂(n) = 0.
OEIS A046109 — Lattice points on a circle of radius n. For the circle of radius n about the origin the entry gives a(n) = 4∏(2eᵢ + 1) over the prime powers pᵢ^eᵢ of n with pᵢ ≡ 1 (mod 4), and a(n) = 8·A046080(n) + 4.
OEIS A046080 — Integer-sided right triangles with hypotenuse n. If fᵢ are the exponents of the primes ≡ 1 (mod 4) in n, then a(n) = (∏(2fᵢ + 1) − 1)/2.
OEIS A004144 — Nonhypotenuse numbers. The sequence 1, 2, 3, 4, 6, 7, 8, 9, 11, 12, … consists, as the entry notes, of the numbers with no prime factor of the form 4k + 1.
Wikipedia — Sieve of Eratosthenes. The multiples of each prime p may be crossed out from p² on; the sieve up to n takes O(n log log n) operations and O(n) memory.
cp-algorithms — Linear Sieve. The linear sieve finds the smallest prime factor of every number up to n in O(n) time and reaches every composite number exactly once.
C++ — Euler233.cpp. The linked immutable C++ revision contains enumerate_core_for_order, tail_fits, build_allowed_multiplier_prefix, solve_sum and run_checkpoints. It is the implementation record for this article.
GCC manual — 128-bit Integers describes the type unsigned __int128. The Java Language Specification, §4.2.1 and §4.2.2: long ranges from −9223372036854775808 to 9223372036854775807, and the integer operators do not indicate overflow or underflow in any way.

Problem 233 source code

C++

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

namespace {

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

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

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

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

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

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

    value = parsed;
    return true;
}

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

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

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

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

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

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

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

    return true;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return 4ULL * product;
}

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

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

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

    return count;
}

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

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

    return spf;
}

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

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

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

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

    return sum;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return true;
}

}  // namespace

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

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

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

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

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

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

    return 0;
}

Python

import math
import multiprocessing
import itertools

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

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

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

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

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

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

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

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

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

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

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

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

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

Java

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return String.valueOf(total);
    }

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