Problem 360: Scary Sphere

Project Euler Problem 360 Solution

EulerSolve provides an optimized solution for Project Euler Problem 360, Scary Sphere, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Define $$S(r)=\sum_{x^2+y^2+z^2=r^2}\left(|x|+|y|+|z|\right),$$ where the sum runs over all integer lattice points on the sphere of radius $r$. For Project Euler 360 the target radius is $r=10^{10}=2^{10}5^{10}$. A direct enumeration of all triples $(x,y,z)$ is far too expensive, so the solution rewrites the three-dimensional problem as a one-dimensional scan combined with the classical sum-of-two-squares function. Mathematical Approach The repository solutions all use the same two ideas: $$\text{(i) symmetry reduces the sum to }r_2(r^2-x^2), \qquad \text{(ii) }10^{10}=2^{10}5^{10}\text{ lets us separate the powers of }2\text{ and }5.$$ Step 1: Reduce the Three-Variable Sum to One Coordinate The three coordinates play identical roles, so the total contribution of $|x|$, $|y|$, and $|z|$ is the same. It is therefore enough to compute the $|x|$-part and multiply by $3$. For a fixed positive $x$, every pair $(y,z)\in\mathbb Z^2$ satisfying $$y^2+z^2=r^2-x^2$$ appears twice in the sphere sum, once with coordinate $x$ and once with coordinate $-x$. Hence the total contribution coming from the $x$-coordinate is $$2x\,r_2(r^2-x^2),$$ where $r_2(n)$ denotes the number of ordered integer pairs $(u,v)$ with $u^2+v^2=n$....

Detailed mathematical approach

Problem Summary

The Manhattan distance between two points of space is the sum of the absolute differences of their coordinates, $|x_1-x_2|+|y_1-y_2|+|z_1-z_2|$. ^[2] Let $C(r)$ be the sphere of radius $r$ centred at the origin $O(0,0,0)$, $I(r)$ the set of all points with integer coordinates on its surface, and $S(r)$ the sum of the Manhattan distances of the elements of $I(r)$ to the origin. The statement gives $S(45)=34518$ and asks for $S(10^{10})$. ^[1]

A point of $I(r)$ is an integer triple with $x^2+y^2+z^2=r^2$ at Manhattan distance $|x|+|y|+|z|$ from the origin, so the quantity wanted is

$$S(r)=\sum_{x^2+y^2+z^2=r^2}\bigl(|x|+|y|+|z|\bigr),$$

the sum running over all such triples.

The whole solution in one sentence: halve the sphere ten times, until its radius is $5^{10}$; then use the symmetry of the three coordinates to add up one coordinate at a time, and count the matching values of the other two from three tables of prime-factor information prepared in advance by a sieve.

Mathematical Approach

1. Four small spheres counted by hand

For $r=1$ the equation $x^2+y^2+z^2=1$ has the six solutions $(\pm1,0,0)$, $(0,\pm1,0)$ and $(0,0,\pm1)$, each at distance $1$, so $S(1)=6$. For $r=2$ only $4=4+0+0$ is possible, six points at distance $2$, so $S(2)=12$. For $r=3$ there are two decompositions of $9$: $9+0+0$ gives $3$ positions for the nonzero coordinate and $2$ signs, six points at distance $3$; $1+4+4$ gives $3$ positions for the $1$ and $2^3=8$ sign patterns, twenty-four points at distance $1+2+2=5$. Hence

$$S(3)=6\cdot3+24\cdot5=18+120=138.$$

For $r=4$ only $16+0+0$ works, so $S(4)=24$.

The start $6,12,138,24$ already shows that $S$ is nowhere near monotone: a radius whose square has many decompositions, such as $3$, dwarfs its larger neighbour $4$. The lattice points really have to be counted, and at $r=10^{10}$ there are $58\,593\,750$ of them.

2. Collapsing the sum to one coordinate

The equation $x^2+y^2+z^2=r^2$ is symmetric in the three coordinates, so exchanging two coordinates is a bijection of $I(r)$ onto itself. The three sums $\sum|x|$, $\sum|y|$ and $\sum|z|$ over $I(r)$ are therefore equal, and

$$S(r)=3\sum_{(x,y,z)\in I(r)}|x|.$$

Now group the points by their first coordinate. For a fixed $x$ with $|x|\le r$, the points of $I(r)$ with that first coordinate are the integer pairs $(y,z)$ with $y^2+z^2=r^2-x^2$. Write

$$r_2(n)=\#\bigl\{(u,v)\in\mathbb Z^2:\ u^2+v^2=n\bigr\},$$

counting ordered pairs, with signs and zeros allowed. The slice $x$ contributes $|x|\,r_2(r^2-x^2)$, the slices $x$ and $-x$ contribute equally, and $x=0$ contributes nothing. Multiplying by $3$ for the coordinates and by $2$ for the sign of $x$ gives the identity the whole solution rests on:

$$S(r)=6\sum_{x=1}^{r}x\,r_2(r^2-x^2).$$

3. The two-squares counting function

For a positive integer $n$, Jacobi's two-square theorem evaluates $r_2(n)$ from its prime factorisation. Write

$$n=2^{c}\prod_i p_i^{e_i}\prod_j q_j^{f_j},$$

where the $p_i$ are the primes congruent to $1$ modulo $4$ and the $q_j$ those congruent to $3$ modulo $4$. Then $r_2(n)=0$ as soon as some $f_j$ is odd, and otherwise

$$r_2(n)=4\prod_i(e_i+1).$$

Two features matter below: the exponent $c$ of $2$ does not appear, and one prime $q\equiv3$ with an odd exponent annihilates the value. ^[3] ^[4]

The reason behind the formula is unique factorisation among the Gaussian integers, the numbers $u+vi$ with $u,v\in\mathbb Z$ and $i^2=-1$. Multiplying $u+vi$ by its conjugate $u-vi$ gives its norm $N(u+vi)=u^2+v^2$, and conjugation respects multiplication, so norms multiply and $r_2(n)$ is the number of Gaussian integers of norm $n$. Factorisation is unique up to order and up to the four numbers $\pm1,\pm i$ of norm $1$, the units, which supply the factor $4$. ^[5] A prime $p\equiv1$ splits into two conjugate factors of norm $p$, as in $5=(2+i)(2-i)$, so a norm $p^{e}$ may take $j$ copies of one and $e-j$ of the other: $e+1$ choices, and the units turn $2\pm i$ into the $r_2(5)=8$ Gaussian integers of norm $5$. A prime $q\equiv3$ stays prime with norm $q^2$, so its exponent in a norm must be even and is then forced, and $2=-i(1+i)^2$ with $N(1+i)=2$ fixes the exponent of $1+i$ with no choice at all, which is why $c$ is absent.

At $n=0$ the product formula does not apply: the only pair is $(0,0)$, so $r_2(0)=1$, the constant term of the theta series of the square lattice. ^[4] Small values used below: $r_2(16)=4$, $r_2(9)=4$, $r_2(25)=12$, and $r_2(24)=r_2(21)=0$ because $24=2^3\cdot3$ and $21=3\cdot7$ each carry a prime $3$ modulo $4$ to an odd power.

4. Worked example: the radius 5

With $r=5$ the five terms of Section 2 can be checked by hand. For each $x$ we factorise $25-x^2$ and read off $r_2$:

$x=1$: $24=2^3\cdot3$, the exponent of $3$ is odd, so $r_2=0$;
$x=2$: $21=3\cdot7$, both exponents odd, so $r_2=0$;
$x=3$: $16=2^4$, no odd prime at all, so $r_2=4$ and the term is $6\cdot3\cdot4=72$;
$x=4$: $9=3^2$, an even exponent and no prime $1$ modulo $4$, so $r_2=4$ and the term is $96$;
$x=5$: $25-25=0$, so $r_2=1$ and the term is $6\cdot5\cdot1=30$.

Hence $S(5)=72+96+30=198$. A direct enumeration of the sphere of radius $5$ finds the same value: its $30$ lattice points are the six signed permutations of $(5,0,0)$, at distance $5$, and the twenty-four signed permutations of $(3,4,0)$, at distance $7$; and $6\cdot5+24\cdot7=198$.

5. Why the formula cannot be used directly

At $r=10^{10}$ the sum of Section 2 has $10^{10}$ terms, each needing the full factorisation of a number up to twenty digits long; trial division alone would cost some $10^{20}$ operations, and the arguments are spread over the whole range up to $10^{20}$, so no table of factorisations can be built either.

The identity $r^2-x^2=(r-x)(r+x)$ shows the way out: the halves never exceed $2r$, a range a sieve can handle. The next two sections make this work, first by shrinking $r$, then by treating the halves independently.

6. Peeling off the powers of two

Every square is $0$ or $1$ modulo $4$, so three of them add up to $0$ modulo $4$ only as $0+0+0$; the other possible sums are $1$, $2$ and $3$. Let $R=2s$ be even. Then $R^2=4s^2\equiv0$, so every point of $I(R)$ has all three coordinates even, and halving them is a bijection from $I(R)$ onto $I(s)$. It divides every Manhattan distance by $2$, hence

$$S(2s)=2\,S(s),\qquad |I(2s)|=|I(s)|.$$

The second identity is an instance of the classical $r_3(4n)=r_3(n)$ for representations as a sum of three squares, which the same argument proves. ^[6] Applying the first ten times to $10^{10}=2^{10}5^{10}$ gives

$$S(10^{10})=2^{10}\,S\bigl(5^{10}\bigr),$$

which replaces the radius $10^{10}$ by $r=5^{10}=9\,765\,625$. This step is what makes the sieve affordable: tables up to $2\cdot10^{10}$ would need about $200$ GB at ten bytes per entry, while tables up to $2\cdot5^{10}=19\,531\,250$ need about $195$ MB.

7. A difference of squares with almost coprime halves

From now on $r=5^{k}$. For $1\le x\lt r$ put

$$a=r-x,\qquad b=r+x,\qquad ab=r^2-x^2.$$

Any common divisor of $a$ and $b$ divides both $a+b=2r$ and $b-a=2x$, hence divides $\gcd(2r,2x)=2\gcd(r,x)$, and $\gcd(r,x)$ divides $r=5^{k}$. Therefore

$$\gcd(a,b)\mid 2\cdot5^{k},$$

so the only primes that can divide both halves are $2$ and $5$. The bound is not far off: for $r=5^4$ the largest common divisor that actually occurs is $250=2\cdot5^3$, at $x=125$.

8. Three tables that give every value of $r_2$

Let $v_p(m)$ denote the exponent of the prime $p$ in $m$. Define, for every positive integer $m$,

$$\begin{aligned}\nu(m)&=v_5(m),\\ g(m)&=\prod_{p\equiv1\,(4),\ p\ne5}\bigl(v_p(m)+1\bigr),\end{aligned}$$

and call $m$ bad when some prime $q\equiv3$ modulo $4$ has an odd exponent in $m$. The claim is that for $1\le x\lt r=5^{k}$, with $a$ and $b$ as in Section 7, the value $r_2(ab)$ is $0$ when $a$ or $b$ is bad and otherwise

$$r_2(ab)=4\,g(a)\,g(b)\bigl(\nu(a)+\nu(b)+1\bigr).$$

The proof uses one tiny observation: if $u$ and $v$ are nonnegative integers with $uv=0$, then

$$(u+1)(v+1)=u+v+1.$$

Exponents add over a product, $v_p(ab)=v_p(a)+v_p(b)$, and by Section 7 at most one of $v_p(a)$, $v_p(b)$ is nonzero for every odd $p\ne5$. Take the three kinds of prime in turn.

A prime $q\equiv3$ modulo $4$ divides at most one half, so $v_q(ab)$ is odd exactly when $v_q(a)$ or $v_q(b)$ is odd: the halves cannot cancel. Hence $ab$ is bad precisely when $a$ or $b$ is, which by Section 3 is the case $r_2(ab)=0$.
A prime $p\equiv1$ modulo $4$ with $p\ne5$ also divides at most one half, so $v_p(a)v_p(b)=0$ and the observation gives $v_p(ab)+1=(v_p(a)+1)(v_p(b)+1)$; the product over all such $p$ is $g(a)g(b)$.
The prime $5$ is congruent to $1$ modulo $4$ and is the one odd prime that may divide both halves, so its factor must come from the sum of the exponents: $v_5(ab)+1=\nu(a)+\nu(b)+1$. Hence its own table.

The prime $2$ needs no treatment, since its exponent does not occur in Jacobi's formula, and multiplying the three contributions proves the claim. Once the tables are built, a slice costs at most six array lookups and a few arithmetic operations — fewer when the bad flags reject it at once — and no further factorisation.

9. Worked example: the radius 25

Take $k=2$, so $r=25$. Of the $24$ values $1\le x\lt25$ only seven survive the bad test; for those, the tables give:

$x=7$: $a=18=2\cdot3^2$, $b=32=2^5$; $g=1,1$ and $\nu=0,0$, so $r_2=4$, term $168$;
$x=9$: $a=16$, $b=34=2\cdot17$; $g(b)=2$, so $r_2=8$, term $432$;
$x=12$: $a=13$, $b=37$, both prime and $1$ modulo $4$; $g=2,2$, so $r_2=16$, term $1152$;
$x=15$: $a=10=2\cdot5$, $b=40=2^3\cdot5$; here $\nu=1,1$, so $r_2=4\cdot3=12$, term $1080$;
$x=16$: $a=9=3^2$, $b=41$; the even exponent of $3$ is harmless and $g(b)=2$, so $r_2=8$, term $768$;
$x=20$: $a=5$, $b=45=3^2\cdot5$; again $\nu=1,1$, so $r_2=12$, term $1440$;
$x=24$: $a=1$, $b=49=7^2$; the even exponent of $7$ is harmless, so $r_2=4$, term $576$.

The endpoint $x=25$ adds $6\cdot25\cdot r_2(0)=150$, and the total is

$$\begin{aligned}S(25)&=168+432+1152+1080\\&\quad+768+1440+576+150\\&=5766.\end{aligned}$$

At $x=15$ and at $x=20$ each half carries one factor $5$. Folding $5$ into $g$ would multiply two factors $1+1$ and give $4$, hence $r_2=16$; since exponents add in the product, the correct factor is $1+1+1=3$ and $r_2=12$. Hence the separate table. Section 6 also gives $S(50)=2\cdot5766=11532$.

10. How far the tables must reach, and how small their entries are

For $1\le x\lt r=5^{k}$ the half $a=r-x$ runs through $1,\dots,r-1$ and the half $b=r+x$ through $r+1,\dots,2r-1$, so indices up to $L=2r$ are enough; the code allocates $0,\dots,L$, that is $L+1=19\,531\,251$ entries at $k=10$. The endpoint $x=r$ stays outside the loop, because $r_2(0)=1$ is not given by Jacobi's formula; it contributes $6r$. The entries are small. Since $5^{11}\gt L$, the largest $\nu(m)$ with $m\le L$ is $10$, attained at $m=5^{10}$ and at $m=2\cdot5^{10}=L$, so one byte holds it. For the other table, $e+1\le p^{e}$ for every prime $p$ and every $e\ge1$; multiplying these inequalities over the primes counted by $g$ gives $g(m)\le m\le L$, well inside a 32-bit word. An exhaustive scan gives the sharper maximum $g(m)=32$, attained for instance at $m=13\cdot17\cdot29\cdot37\cdot41$, that is $9\,722\,453$.

11. The finished recipe

Collecting Sections 2, 6, 8 and 10, the computation is

$$S(10^{10})=2^{10}\Bigl(6r+6\sum_{x=1}^{r-1}x\,r_2(ab)\Bigr),$$

with $r=5^{10}$, $a=r-x$, $b=r+x$ and $r_2(ab)$ from Section 8; the rejected slices need no exclusion from the sum, since their $r_2(ab)$ is $0$. Only $475\,088$ of the $9\,765\,624$ values of $x$ have both halves good, about $4.86$ per cent of them, so most iterations do nothing but two table lookups and a branch. Among the surviving terms the largest value of $r_2(r^2-x^2)$ is $1792$, the largest $g$ that occurs is $32$ and the largest $\nu(a)+\nu(b)$ is $18$.

How the Code Works

The C++ program is organised like the sections above. pow_u64 builds $r=5^{k}$ by repeated squaring, and isqrt_u64 is an integer square root used only by the brute-force check. build_spf is a sieve of Eratosthenes that, instead of a Boolean flag, stores in spf[n] the smallest prime factor of $n$: it initialises spf[n] = n and, for every $i$ that is still its own smallest factor, overwrites the untouched multiples $i^2, i^2+i,\dots$ with $i$. Any $n$ can then be factorised by repeatedly dividing by spf[n].

build_factor_tables walks $n=1,\dots,L$ and performs that factorisation, filling the three tables of Section 8 in one pass: bad_3mod4_odd[n] is set when a prime $p\equiv3$ modulo $4$ comes out with an odd exponent, exp5[n] receives the exponent of $5$, and one_mod4_factor[n] accumulates the product of $e+1$ over the primes $p\equiv1$ modulo $4$ with $p\ne5$; the two residue classes are tested as (p & 3) == 3 and (p & 3) == 1. The narrow types u8 and u32 are justified by Section 10.

sum_for_radius_power_of_five builds the tables for $L=2r$ and then runs the loop of Section 11 over $x=1,\dots,r-1$, skipping $x$ as soon as bad[a] or bad[b] is set and otherwise adding $6x\cdot$ways with ways $=4\,g(a)g(b)(\nu(a)+\nu(b)+1)$. The range of $x$ is cut into contiguous blocks, each keeping its own 128-bit partial sum; the partial sums are then added to the endpoint term $6r$, the $x=r$ contribution $r_2(0)=1$. Integer addition is exact and associative, so the result does not depend on how the range is cut. sum_for_radius_2a_5b adds the last step of Section 6 as the shift total <<= exp2, and main calls it with $(10,10)$, printing the value through to_string_u128, since the standard streams do not print the 128-bit type.

The accumulators are deliberately wide. Section 2 on its own bounds the sum only crudely: $r_2(n)\le4d(n)$ with $d$ the number of divisors, and $d(n)\le2\sqrt n$ since divisors pair as $d$ and $n/d$, so each term is at most $48r^2$ and the sum at most $48r^3$, about $5\cdot10^{31}$ at $r=10^{10}$ — far outside the $64$-bit range. C++ therefore converts the table entries to unsigned __int128 before multiplying them into ways. ^[9]

brute_sum is the independent check: it loops over $x$ and $y$ in $[-r,r]$, tests whether $r^2-x^2-y^2$ is a perfect square, and adds $|x|+|y|$ once when $z=0$ and $|x|+|y|+|z|$ twice otherwise, the correct treatment of the two signs of $z$. It computes $r^2$ in int, so it suits only the small radii it is called with. run_checkpoints compares brute_sum(45) with $34518$ and the fast route with $S(5)=198$, $S(10)=396$ and $S(25)=5766$; main runs it before the target computation and aborts with exit code $2$ on failure, unless given --skip-checkpoints.

The two ports do the same mathematics in a different layout. Python puts the sieve, the tables, the summation and the final scaling inside one function solve, keeps the sieve in a list built as list(range(limit + 1)) and relies on unlimited integer precision, so no accumulator width is needed. Java keeps the stages apart in buildSpf, buildFactorTables, sumForRadiusPowerOfFive and solve, holds the tables in int[], boolean[] and byte[], computes waysVal in a long — safe because $4g(a)g(b)(\nu(a)+\nu(b)+1)$ stays below a few thousand — and accumulates in BigInteger, finishing with shiftLeft(10) for the factor $2^{10}$. ^[9] It obtains $r$ as (long) Math.pow(5, exp), through a floating-point power; the value is exact here, since a double represents $5^{10}$ exactly, but an integer loop would not depend on the accuracy of Math.pow. Neither port runs the checkpoints; only the C++ program does.

Complexity and Verification

Let $L=2r=19\,531\,250$ for the target. The sieve costs $O(L\log\log L)$ time and $O(L)$ memory. ^[7] Filling the tables divides each $n\le L$ down to $1$, which takes $\Omega(n)$ divisions, with $\Omega$ counting prime factors with multiplicity; over the range that is

$$\sum_{n\le L}\Omega(n)=L\log\log L+B_2L+o(L),$$

with $B_2=1.0345061\ldots$. ^[8] Building the table for the target alone performs exactly $74\,766\,832$ such divisions, $3.8281$ per number, against the predicted $\log\log L+B_2=3.8551$; its outer loop over distinct primes runs $59\,666\,961$ times, $3.0549$ per number, against the predicted $\log\log L+B_1=3.0821$ with the Mertens constant $B_1$ (the default run also builds the small checkpoint tables). The main loop is $O(r)$ with constant work per step, so the whole method is $O(r\log\log r)$ time and $O(r)$ memory, against $O(r^2)$ for brute_sum.

The four tables occupy $4+4+1+1=10$ bytes per index, about $195$ MB for $L+1$ entries. On the machine used for this article the C++ program peaks at about $198$ MB and finishes in under a second, as does the Java port, while Python takes about ten seconds and about $1.27$ GB; both figures depend on the machine, the compiler or interpreter and the measurement method. Python holds four lists of object references instead of packed arrays of machine words, and list(range(limit + 1)) first creates $L+1$ separate integer objects.

Beyond the packaged checkpoints, the following were verified for this article:

A triple loop over the whole sphere and the formula $6\sum x\,r_2(r^2-x^2)$ of Section 2 agree for every $1\le r\le30$ and for $r=32,40,45,50,55,60$${},64,65,75,80,100,125$.
Jacobi's formula for $r_2(n)$ agrees with a direct count of the pairs $(u,v)$ for every $0\le n\lt3000$, including $r_2(0)=1$.
The product identity of Section 8 agrees with a direct factorisation of $r^2-x^2$ for every $1\le x\lt5^4$, and every $\gcd(r-x,r+x)$ there is of the form $2^i5^j$ dividing $2\cdot5^4$, as Section 7 requires.
$S(2s)=2S(s)$ holds for every $1\le s\le60$, and the $64$ residue triples modulo $4$ confirm the descent of Section 6.
The table method reproduces $S(5)=198$, $S(25)=5766$, $S(125)=142038$, $S(625)=3522390$ and $S(3125)=87758214$, each also obtained directly from the factorisation formula.
The tables give $|I(10^{10})|=|I(5^{10})|=$${}58\,593\,750=6\cdot5^{10}$ lattice points, with $r_2(5^{20})=4\cdot21=84$ of them in the slice $x=0$. A direct count confirms $|I(5^{k})|=6\cdot5^{k}$ for $k\le6$ and $|I(45)|=510$; that pattern is an observation, not a step of the solution, but it cross-checks the tables at full size.
The three published programs were run to completion at the target radius and printed the same value.

Footnotes and References

The references below supply the background facts used above. The reduction to one coordinate, the descent that removes the powers of two, the near-coprimality of $r-x$ and $r+x$ and the splitting of $r_2$ into three tables are derived explicitly in this article.

Project Euler 360 — Scary Sphere. The official statement defines the Manhattan distance between two points as |x₁ − x₂| + |y₁ − y₂| + |z₁ − z₂|, lets C(r) be the sphere of radius r centred at the origin, I(r) the set of its surface points with integer coordinates and S(r) the sum of their Manhattan distances to the origin, gives S(45) = 34518 and asks for S(10¹⁰).
Wikipedia — Taxicab geometry. It defines the taxicab or Manhattan distance as the sum of the absolute differences of the respective Cartesian coordinates, gives the n-dimensional formula as the sum of |pᵢ − qᵢ| and identifies it with the L¹ distance.
Wikipedia — Sum of squares function. It defines r_k(n) as the number of ordered k-tuples of integers whose squares sum to n, with signs counted separately and zero allowed, and states the factorisation form of Section 3: for n = 2^g · ∏pᵢ^fᵢ · ∏qⱼ^hⱼ with pᵢ ≡ 1 and qⱼ ≡ 3 modulo 4, r₂(n) is 4∏(fᵢ + 1) if every hⱼ is even and 0 otherwise. It also gives the equivalent divisor form 4(d₁(n) − d₃(n)).
OEIS A004018 — theta series of the square lattice. The entry lists the values of r₂, beginning 1, 4, 4, 0, 4, 8, 0, 0, 4, 4, 8 at n = 0, 1, 2, …, states a(0) = 1 and repeats the factorisation formula of Section 3. Its first terms confirm r₂(0) = 1, r₂(9) = 4, r₂(16) = 4 and r₂(21) = r₂(24) = 0.
Wikipedia — Gaussian integer. It states that the norm N(a + bi) = a² + b² is multiplicative, that ℤ[i] is a unique factorisation domain whose units are 1, −1, i, −i, and that a rational prime p ≡ 3 modulo 4 stays prime while p ≡ 1 modulo 4 splits, as in 5 = (2 + i)(2 − i), with 2 = (1 + i)(1 − i) ramified — the facts used in Section 3.
OEIS A005875 — theta series of the simple cubic lattice. The entry counts the ordered integer triples whose squares sum to n, beginning 1, 6, 12, 8, 6, 24 at n = 0, 1, …, and records a(4n) = a(n) and a(8n + 7) = 0. Its a(25) = 30 is the point count used in Section 4, and a(4n) = a(n) is the counting half of the descent of Section 6.
Wikipedia — Sieve of Eratosthenes. It gives the sieve’s time as O(n log log n) operations and its memory as O(n). The page describes the plain Boolean sieve; the variant used here, which records the smallest prime factor instead of a flag and so doubles as a factorisation table, is derived in the article itself.
Wikipedia — Prime omega function. It defines ω(n) as the number of distinct prime factors and Ω(n) as the number with multiplicity, and gives the average orders used above: the sum of ω(n) over n ≤ x is x log log x + B₁x + o(x) with the Mertens constant B₁ = 0.2614972128…, and the sum of Ω(n) is x log log x + B₂x + o(x) with B₂ = 1.0345061758….
GCC manual — 128-bit integers. It documents unsigned __int128 on targets with an integer mode wide enough to hold 128 bits; the C++ file uses it for the accumulator and in to_string_u128. Java — BigInteger documents the arbitrary-precision integers of the Java port, including the shiftLeft(n) that multiplies by 2ⁿ, and Python — numeric types states that Python integers have unlimited precision.

Problem 360 source code

C++

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

namespace {

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

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

u64 pow_u64(u64 base, unsigned exp) {
    u64 result = 1;
    while (exp > 0U) {
        if ((exp & 1U) != 0U) {
            result *= base;
        }
        base *= base;
        exp >>= 1U;
    }
    return result;
}

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

std::vector<u32> build_spf(const u32 limit) {
    std::vector<u32> spf(static_cast<std::size_t>(limit + 1U));
    for (u32 i = 0; i <= limit; ++i) {
        spf[static_cast<std::size_t>(i)] = i;
    }

    for (u32 i = 2; static_cast<u64>(i) * i <= limit; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (u32 j = i * i; j <= limit; j += i) {
            if (spf[static_cast<std::size_t>(j)] == j) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }

    return spf;
}

struct FactorTables {
    std::vector<u8> bad_3mod4_odd;
    std::vector<u8> exp5;
    std::vector<u32> one_mod4_factor;
};

FactorTables build_factor_tables(const u32 limit, const std::vector<u32>& spf) {
    FactorTables data;
    data.bad_3mod4_odd.assign(static_cast<std::size_t>(limit + 1U), 0U);
    data.exp5.assign(static_cast<std::size_t>(limit + 1U), 0U);
    data.one_mod4_factor.assign(static_cast<std::size_t>(limit + 1U), 1U);

    for (u32 n = 1; n <= limit; ++n) {
        u32 x = n;
        bool bad = false;
        u8 e5 = 0;
        u32 mult = 1;

        while (x > 1U) {
            const u32 p = spf[static_cast<std::size_t>(x)];
            u32 e = 0;
            do {
                x /= p;
                ++e;
            } while (x % p == 0U);

            if (p == 5U) {
                e5 = static_cast<u8>(e);
            } else if ((p & 3U) == 3U) {
                if ((e & 1U) != 0U) {
                    bad = true;
                }
            } else if ((p & 3U) == 1U) {
                mult *= (e + 1U);
            }
        }

        data.bad_3mod4_odd[static_cast<std::size_t>(n)] = bad ? 1U : 0U;
        data.exp5[static_cast<std::size_t>(n)] = e5;
        data.one_mod4_factor[static_cast<std::size_t>(n)] = mult;
    }

    return data;
}

u128 sum_for_radius_power_of_five(const unsigned exponent) {
    const u64 r = pow_u64(5ULL, exponent);
    const u32 limit = static_cast<u32>(2ULL * r);

    const std::vector<u32> spf = build_spf(limit);
    const FactorTables tables = build_factor_tables(limit, spf);

    const auto& bad = tables.bad_3mod4_odd;
    const auto& e5 = tables.exp5;
    const auto& f1 = tables.one_mod4_factor;

    const u64 work_items = (r >= 1ULL) ? (r - 1ULL) : 0ULL;
    unsigned thread_count = 1U;
    if (work_items >= 1000000ULL) {
        thread_count = std::thread::hardware_concurrency();
        if (thread_count == 0U) {
            thread_count = 4U;
        }
        thread_count = static_cast<unsigned>(std::min<u64>(thread_count, work_items));
    }

    std::vector<u128> partial(thread_count, 0U);
    std::vector<std::thread> workers;
    workers.reserve(thread_count);

    const u64 chunk = (work_items + thread_count - 1ULL) / thread_count;
    for (unsigned t = 0; t < thread_count; ++t) {
        const u64 begin = 1ULL + static_cast<u64>(t) * chunk;
        const u64 end = std::min<u64>(r, begin + chunk);

        workers.emplace_back([&, t, begin, end]() {
            u128 local = 0;
            for (u64 x = begin; x < end; ++x) {
                const u32 a = static_cast<u32>(r - x);
                const u32 b = static_cast<u32>(r + x);

                if (bad[static_cast<std::size_t>(a)] != 0U || bad[static_cast<std::size_t>(b)] != 0U) {
                    continue;
                }

                const u128 ways = 4U * static_cast<u128>(f1[static_cast<std::size_t>(a)]) *
                                  static_cast<u128>(f1[static_cast<std::size_t>(b)]) *
                                  static_cast<u128>(static_cast<u32>(e5[static_cast<std::size_t>(a)]) +
                                                    static_cast<u32>(e5[static_cast<std::size_t>(b)]) + 1U);

                local += 6U * static_cast<u128>(x) * ways;
            }
            partial[t] = local;
        });
    }

    for (auto& th : workers) {
        th.join();
    }

    // x = r contributes r_2(0) = 1.
    u128 total = 6U * static_cast<u128>(r);
    for (const u128 v : partial) {
        total += v;
    }
    return total;
}

u128 sum_for_radius_2a_5b(const unsigned exp2, const unsigned exp5) {
    u128 total = sum_for_radius_power_of_five(exp5);
    total <<= exp2;
    return total;
}

u64 brute_sum(const int r) {
    const int rr = r * r;
    u128 total = 0;

    for (int x = -r; x <= r; ++x) {
        for (int y = -r; y <= r; ++y) {
            const int z2 = rr - x * x - y * y;
            if (z2 < 0) {
                continue;
            }

            const int z = static_cast<int>(isqrt_u64(static_cast<u64>(z2)));
            if (z * z != z2) {
                continue;
            }

            const u64 base = static_cast<u64>(std::abs(x)) + static_cast<u64>(std::abs(y));
            if (z == 0) {
                total += base;
            } else {
                total += base + static_cast<u64>(z);
                total += base + static_cast<u64>(z);
            }
        }
    }

    return static_cast<u64>(total);
}

bool run_checkpoints() {
    if (brute_sum(45) != 34518ULL) {
        std::cerr << "Checkpoint failed: S(45)\n";
        return false;
    }

    if (sum_for_radius_2a_5b(0U, 1U) != 198U) {
        std::cerr << "Checkpoint failed: S(5)\n";
        return false;
    }

    if (sum_for_radius_2a_5b(1U, 1U) != 396U) {
        std::cerr << "Checkpoint failed: S(10)\n";
        return false;
    }

    if (sum_for_radius_2a_5b(0U, 2U) != 5766U) {
        std::cerr << "Checkpoint failed: S(25)\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    // 10^10 = 2^10 * 5^10
    const u128 answer = sum_for_radius_2a_5b(10U, 10U);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

import math

def solve():
    # 10^10 = 2^10 * 5^10, so we compute for radius = 5^10 then shift by 2^10
    exp2 = 10
    exp5 = 10
    r = 5 ** exp5
    limit = 2 * r

    # Build SPF sieve
    spf = list(range(limit + 1))
    i = 2
    while i * i <= limit:
        if spf[i] == i:
            for j in range(i*i, limit+1, i):
                if spf[j] == j:
                    spf[j] = i
        i += 1

    bad = [0] * (limit + 1)
    e5 = [0] * (limit + 1)
    f1 = [1] * (limit + 1)

    for n in range(1, limit + 1):
        x = n
        is_bad = False
        exp_5 = 0
        mult = 1
        while x > 1:
            p = spf[x]
            e = 0
            while x % p == 0:
                x //= p
                e += 1
            if p == 5:
                exp_5 = e
            elif p % 4 == 3:
                if e % 2 != 0:
                    is_bad = True
            elif p % 4 == 1:
                mult *= (e + 1)
        bad[n] = 1 if is_bad else 0
        e5[n] = exp_5
        f1[n] = mult

    total = 0
    for x in range(1, r):
        a = r - x
        b = r + x
        if bad[a] or bad[b]:
            continue
        ways = 4 * f1[a] * f1[b] * (e5[a] + e5[b] + 1)
        total += 6 * x * ways

    # x = r: contributes r_2(0)=1 => 6*r
    total += 6 * r
    total <<= exp2

    return str(total)

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

Java

import java.math.BigInteger;

public class Euler360 {

    static int[] buildSpf(int limit) {
        int[] spf = new int[limit + 1];
        for (int i = 0; i <= limit; i++)
            spf[i] = i;
        for (int i = 2; (long) i * i <= limit; i++) {
            if (spf[i] == i) {
                for (int j = i * i; j <= limit; j += i) {
                    if (spf[j] == j)
                        spf[j] = i;
                }
            }
        }
        return spf;
    }

    static class FactorTables {
        boolean[] bad;
        byte[] e5;
        int[] f1;
    }

    static FactorTables buildFactorTables(int limit, int[] spf) {
        FactorTables t = new FactorTables();
        t.bad = new boolean[limit + 1];
        t.e5 = new byte[limit + 1];
        t.f1 = new int[limit + 1];

        for (int n = 1; n <= limit; n++) {
            int x = n;
            boolean bad = false;
            byte exp5 = 0;
            int mult = 1;

            while (x > 1) {
                int p = spf[x];
                int e = 0;
                do {
                    x /= p;
                    e++;
                } while (x % p == 0);

                if (p == 5) {
                    exp5 = (byte) e;
                } else if ((p & 3) == 3) {
                    if ((e & 1) != 0)
                        bad = true;
                } else if ((p & 3) == 1) {
                    mult *= (e + 1);
                }
            }
            t.bad[n] = bad;
            t.e5[n] = exp5;
            t.f1[n] = mult;
        }
        return t;
    }

    static BigInteger sumForRadiusPowerOfFive(int exp) {
        long r = (long) Math.pow(5, exp);
        int limit = (int) (2 * r);

        int[] spf = buildSpf(limit);
        FactorTables tables = buildFactorTables(limit, spf);

        boolean[] bad = tables.bad;
        byte[] e5 = tables.e5;
        int[] f1 = tables.f1;

        BigInteger total = BigInteger.valueOf(r).multiply(BigInteger.valueOf(6));
        BigInteger mult6 = BigInteger.valueOf(6);
        BigInteger mult4 = BigInteger.valueOf(4);

        for (long x = 1; x < r; x++) {
            int a = (int) (r - x);
            int b = (int) (r + x);

            if (bad[a] || bad[b])
                continue;

            long waysVal = 4L * f1[a] * f1[b] * (e5[a] + e5[b] + 1L);
            BigInteger ways = BigInteger.valueOf(waysVal);
            BigInteger add = BigInteger.valueOf(x).multiply(mult6).multiply(ways);
            total = total.add(add);
        }
        return total;
    }

    static BigInteger solve() {
        return sumForRadiusPowerOfFive(10).shiftLeft(10);
    }

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