Detailed mathematical approach

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\). Multiplying by \(3\) for the three symmetric coordinates gives

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

The term \(x=0\) does not matter because it would be multiplied by \(x\). The endpoint \(x=r\) is special: then \(r^2-x^2=0\), so \(r_2(0)=1\) from the single pair \((0,0)\).

Step 2: Use the Sum-of-Two-Squares Function

For each \(x\), write

$r^2-x^2=(r-x)(r+x)=ab,$

with \(a=r-x\) and \(b=r+x\). The problem is now to evaluate \(r_2(ab)\).

If

$n=2^\alpha \prod_{p\equiv 1 \pmod{4}} p^{e_p}\prod_{q\equiv 3 \pmod{4}} q^{f_q},$

then the classical formula for the ordered signed representation count is

$r_2(n)= \begin{cases} 0, & \text{if some } f_q \text{ is odd},\\ 4\displaystyle\prod_{p\equiv 1 \pmod{4}}(e_p+1), & \text{otherwise}. \end{cases}$

So primes congruent to \(3 \pmod{4}\) are only a parity test: any odd exponent kills the term. Primes congruent to \(1 \pmod{4}\) contribute multiplicative factors \(e_p+1\). The exponent of \(2\) does not affect the value of \(r_2(n)\), so the implementation does not need a separate table for powers of \(2\).

Step 3: Why the Factor Tables for \(r=5^b\) Are Enough

The code first handles radii of the form \(r=5^b\). For \(x\in\{1,\dots,r-1\}\), set

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

Any common divisor of \(a\) and \(b\) divides both \(a+b=2r\) and \(b-a=2x\), hence

$\gcd(a,b)\mid 2r=2\cdot 5^b.$

This simple fact is the key simplification. It implies that any odd prime \(p\neq 5\) can divide at most one of \(a\) and \(b\). Therefore:

\(1.\) A prime \(q\equiv 3 \pmod{4}\) has odd exponent in \(ab\) exactly when it has odd exponent in \(a\) or in \(b\). It is enough to know whether each factor is “bad”.

\(2.\) A prime \(p\equiv 1 \pmod{4}\), \(p\neq 5\), contributes independently from \(a\) and \(b\), because it cannot appear in both factors.

\(3.\) The prime \(5\) must be treated separately, because it is the only prime congruent to \(1 \pmod{4}\) that may appear in both \(a\) and \(b\).

That is why the solutions precompute, for every \(m\le 2r\):

$\texttt{bad}(m)= \begin{cases} 1,& \text{if some } q\equiv 3 \pmod{4} \text{ has odd exponent in }m,\\ 0,& \text{otherwise}, \end{cases}$

$\nu_5(m)=v_5(m),$

$g(m)=\prod_{\substack{p\equiv 1 \pmod{4}\\ p\neq 5}}(e_p(m)+1).$

If either \(a\) or \(b\) is bad, then \(r_2(ab)=0\). Otherwise the exponent of \(5\) in \(ab\) is \(\nu_5(a)+\nu_5(b)\), and the remaining \(1 \pmod{4}\) primes factor cleanly, so

$\boxed{r_2(ab)=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).}$

This is exactly the expression used in the C++, Python, and Java implementations.

Step 4: Recover the Radius \(2^a5^b\)

Now let the radius be \(R=2^a5^b\). Because \(R^2\) is divisible by \(4\), any solution of

$x^2+y^2+z^2=R^2$

must have all three coordinates even: modulo \(4\), each square is \(0\) or \(1\), and the only way three such residues can sum to \(0\) is \(0+0+0\). Dividing by \(2\) and repeating this argument \(a\) times shows that every solution has the form

$x=2^a x',\qquad y=2^a y',\qquad z=2^a z',$

with

$x'^2+y'^2+z'^2=5^{2b}.$

This is a bijection between lattice points on the sphere of radius \(2^a5^b\) and those on the sphere of radius \(5^b\). The absolute-value sum scales by the same factor \(2^a\), so

$\boxed{S(2^a5^b)=2^aS(5^b).}$

That explains the final left shift in the code after the \(5^b\) case has been computed.

Step 5: Worked Example with \(r=5\)

The small checkpoint \(r=5\) is useful because the whole derivation can be checked by hand:

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

Evaluate each term:

\(x=1\): \(25-1=24=2^3\cdot 3\), so a prime \(3 \pmod{4}\) has odd exponent and \(r_2(24)=0\).

\(x=2\): \(25-4=21=3\cdot 7\), again impossible, so \(r_2(21)=0\).

\(x=3\): \(25-9=16\), so \(r_2(16)=4\) from \((\pm4,0)\) and \((0,\pm4)\).

\(x=4\): \(25-16=9\), so \(r_2(9)=4\) from \((\pm3,0)\) and \((0,\pm3)\).

\(x=5\): \(25-25=0\), so \(r_2(0)=1\).

Therefore

$S(5)=6\cdot 3\cdot 4 + 6\cdot 4\cdot 4 + 6\cdot 5 = 72+96+30=198.$

The scaling law then gives

$S(10)=2S(5)=396,$

which matches the checkpoint used by the implementation.

How the Code Works

The C++ solution builds an SPF table up to \(2r\), where \(r=5^{10}\), with build_spf. Then build_factor_tables factors every integer \(1\le m\le 2r\) and stores the three quantities described above: the bad flag, the exponent of \(5\), and the multiplicative factor from primes congruent to \(1 \pmod{4}\) other than \(5\).

The main summation loops over \(x=1,\dots,r-1\), computes \(a=r-x\) and \(b=r+x\), rejects the pair immediately if either factor is bad, and otherwise evaluates

$\texttt{ways}=4\,g(a)g(b)\bigl(\nu_5(a)+\nu_5(b)+1\bigr).$

It then adds

$6x\cdot \texttt{ways}$

to the total. After the loop it adds the endpoint contribution \(6r\) for \(x=r\), and finally multiplies by \(2^{10}\) via a left shift.

The Python and Java files implement the same arithmetic formula. The Java version stores the final answer in BigInteger. The C++ version uses unsigned __int128 and parallelizes the main \(x\)-range across worker threads. It also contains explicit checkpoints: \(S(5)=198\), \(S(10)=396\), \(S(25)=5766\), and a brute-force verification \(S(45)=34518\).

Complexity Analysis

Let \(r=5^{10}\) and \(\texttt{limit}=2r\). Building the SPF sieve costs \(O(\texttt{limit}\log\log \texttt{limit})\) time and \(O(\texttt{limit})\) memory. Filling the factor tables by repeatedly stripping SPF factors is near-linear on average, again about \(O(\texttt{limit}\log\log \texttt{limit})\) time. The final scan over \(x=1,\dots,r-1\) is \(O(r)\) with only constant-time table lookups per iteration. The overall memory usage is \(O(r)\), dominated by the precomputed arrays. The C++ threading changes wall-clock time, but not the asymptotic complexity.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=360
2. Sum of two squares theorem: Wikipedia — Sum of two squares theorem
3. Sum of two squares function: Wikipedia — Sum of two squares function
4. Sieve of Eratosthenes and SPF tables: Wikipedia — Sieve of Eratosthenes