Problem 273: Sum of Squares

Detailed mathematical approach

Problem Summary

For a positive integer \(N\), consider all pairs of integers \((a,b)\) with \(a^2+b^2=N\) and \(0\le a\le b\), and let \(S(N)\) be the sum of the values of \(a\) over these pairs. The statement gives \(65=1^2+8^2=4^2+7^2\), so that \(S(65)=1+4=5\). We must find the sum of \(S(N)\) over all squarefree \(N\), that is, \(N\) not divisible by the square of any prime, that are only divisible by primes of the form \(4k+1\) below \(150\). ^[1] There are \(16\) such primes, 5, 13, 17, 29, 37, 41, 53, 61, 73, 89, 97, 101, 109, 113, 137 and 149, so the sum runs over the \(2^{16}-1=65535\) products of one or more distinct primes from this list; the empty product \(N=1=0^2+1^2\) would only add \(S(1)=0\).

The whole solution in one sentence: every prime of the form \(4k+1\) is the product of two conjugate complex numbers with integer parts, \(u+v\,i\) and \(u-v\,i\); multiplying one of the two factors of every prime of \(N\), with the factor of the smallest prime held fixed, produces each solution of \(a^2+b^2=N\) exactly once, so the program adds up the smaller part, without sign, of all these products, about \(21.5\) million of them.

Mathematical Approach

We start from small cases, produce all solutions by multiplying Gaussian integers, prove that each solution appears exactly once and check that 64-bit integers suffice.

1. Small cases

A prime of the list has one solution: \(5=1^2+2^2\), \(13=2^2+3^2\) and \(17=1^2+4^2\), so \(S(5)=1\), \(S(13)=2\) and \(S(17)=1\). A product of two of these primes has two solutions, and the product of all three has four:

• \(65=1^2+8^2=4^2+7^2\), so \(S(65)=5\);
• \(85=2^2+9^2=6^2+7^2\), so \(S(85)=8\);
• \(221=5^2+14^2=10^2+11^2\), so \(S(221)=15\);
• \(1105=4^2+33^2=9^2+32^2\)\({}=12^2+31^2=23^2+24^2\), so \(S(1105)=48\).

If the list consisted only of 5, 13 and 17, the requested sum would be \(1+2+1+5+8+15+48\)\({}=80\). The number of solutions doubles with every new prime factor. Moreover, \(a=0\) would make \(N=b^2\) a square, which a squarefree \(N\gt1\) is not, and \(a=b\) would make \(N=2a^2\) even, while our primes are odd; so every solution has \(0\lt a\lt b\). A direct search is hopeless: for the largest \(N\), the product of all \(16\) primes, about \(2.5\cdot10^{27}\), the value of \(a\) would have to run up to about \(3.5\cdot10^{13}\).

2. Which primes are sums of two squares

Every square leaves the remainder \(0\) or \(1\) on division by \(4\), since \((2m)^2=4m^2\) and \((2m+1)^2=4m^2+4m+1\), so a sum of two squares never has the form \(4k+3\). Fermat's theorem on sums of two squares says that every prime \(p\equiv1\pmod4\), a so-called Pythagorean prime, is a sum of two squares. ^{[2] [3]} As \(p\) is odd, exactly one of the two squares is odd, and we write \(p=u_p^2+v_p^2\) with \(u_p\) odd and \(v_p\) even, both positive:

\(5=1^2+2^2\), \(13=3^2+2^2\), \(17=1^2+4^2\), \(29=5^2+2^2\), \(37=1^2+6^2\), \(41=5^2+4^2\), \(53=7^2+2^2\), \(61=5^2+6^2\), \(73=3^2+8^2\), \(89=5^2+8^2\), \(97=9^2+4^2\), \(101=1^2+10^2\), \(109=3^2+10^2\), \(113=7^2+8^2\), \(137=11^2+4^2\) and \(149=7^2+10^2\).

Section 5 shows that these representations are unique.

3. Gaussian integers and the product rule

A Gaussian integer is a complex number \(z=x+y\,i\) whose parts \(x\) and \(y\) are integers, positive, negative or zero, where \(i^2=-1\). ^[4] Its conjugate is \(\overline z=x-y\,i\), and its norm is \(|z|^2=z\overline z=x^2+y^2\). So every Gaussian integer of norm \(N\) gives a solution of \(a^2+b^2=N\): its two parts, taken without signs and sorted, are \(a\) and \(b\). Products are computed with the usual rule

$$\begin{aligned}(x+y\,i)(u+v\,i)&=(xu-yv)\\&\quad+(xv+yu)\,i.\end{aligned}$$

Since \(|zw|=|z|\,|w|\) for complex numbers, the norm is multiplicative. Written out, this is the Brahmagupta–Fibonacci identity ^[5]

$$\begin{aligned}(x^2+y^2)(u^2+v^2)&=(xu-yv)^2\\&\quad+(xv+yu)^2,\end{aligned}$$

and multiplying by the conjugate, which replaces \(v\) by \(-v\), gives the second representation \((xu+yv)^2+(xv-yu)^2\). We write \(z_p=u_p+v_p\,i\), so that \(p=z_p\overline{z_p}\). For the example of the statement, \(z_5=1+2i\) and \(z_{13}=3+2i\) give

$$\begin{aligned}(1+2i)(3+2i)&=-1+8i,\\(1+2i)(3-2i)&=7+4i,\end{aligned}$$

exactly the two solutions \(65=1^2+8^2=4^2+7^2\).

4. Eight Gaussian integers for one solution

The units, the Gaussian integers whose inverse is again a Gaussian integer, are the numbers of norm \(1\): \(\pm1\) and \(\pm i\). ^[4] Multiplying \(x+y\,i\) by \(i\) gives \(-y+x\,i\): the two parts swap places and one sign changes. So the four associates of \(z\), the numbers \(z\), \(-z\), \(iz\) and \(-iz\), describe the same solution, and so do the four associates of \(\overline z\). For a solution with \(0\lt a\lt b\) these are the eight different numbers \(\pm a\pm b\,i\) and \(\pm b\pm a\,i\), and no other Gaussian integer describes it.

5. Gaussian primes and unique factorisation

A Gaussian prime is a Gaussian integer, other than \(0\) and the units, that is not the product of two non-units. Gaussian integers can be divided with a remainder of smaller norm: to divide \(\alpha\) by \(\beta\ne0\), rounding both parts of \(\alpha/\beta\) to the nearest integers gives a Gaussian integer \(q\) with \(|\alpha/\beta-q|^2\le1/2\), so \(r=\alpha-q\beta\) has \(|r|^2\le|\beta|^2/2\). For example, \(7=(1+2i)(1-3i)+i\), where \(1-3i\) is \(7/(1+2i)=1.4-2.8i\) rounded. As for ordinary integers, this division leads to unique factorisation: every Gaussian integer other than \(0\) and the units is a product of Gaussian primes, unique up to the order of the factors and up to units. ^[4] In particular, a Gaussian prime that divides a product \(\alpha\beta\) divides \(\alpha\) or \(\beta\), because up to a unit it must occur among the prime factors of \(\alpha\) and \(\beta\). Three facts follow.

• \(z_p\) and \(\overline{z_p}\) are Gaussian primes: if \(z_p=\alpha\beta\), then \(|\alpha|^2|\beta|^2=p\), so one factor has norm \(1\) and is a unit.
• \(z_p\) and \(\overline{z_p}\) are not associates: the associates of \(u+v\,i\) are \(u+v\,i\), \(-u-v\,i\), \(-v+u\,i\) and \(v-u\,i\), and \(u-v\,i\) is one of them only if \(v=0\), \(u=0\), \(u=-v\) or \(u=v\). Then \(p=u^2+v^2\) would be a square or twice a square, which an odd prime is not.
• The representation of Section 2 is unique: if \(p=s^2+t^2\), then \(s+t\,i\) is a Gaussian prime, as in the first point, dividing \(p=z_p\overline{z_p}\), hence an associate of \(z_p\) or \(\overline{z_p}\), so \(\{|s|,|t|\}=\{u_p,v_p\}\) by Section 4.

6. Every solution twice, then once

Let \(N=p_1p_2\cdots p_k\) with distinct primes from the list and \(p_1\) the smallest of them. Choosing one factor \(\eta_j\in\{z_{p_j},\overline{z_{p_j}}\}\) for every prime and multiplying gives a Gaussian integer of norm \(N\) by Section 3, hence a solution. There are \(2^k\) such choices, but they come in pairs: switching every choice to the other factor turns the product \(Z\) into \(\overline Z\), which describes the same solution. Exactly one product of each pair uses \(z_{p_1}\), so we fix that factor and keep the \(2^{k-1}\) products

$$Z=z_{p_1}\,\eta_2\,\eta_3\cdots\eta_k,$$

in which every \(\eta_j\) with \(j\ge2\) is still either \(z_{p_j}\) or \(\overline{z_{p_j}}\). For \(N=65\) these are the two products shown in Section 3.

7. The key claim and its proof

Claim. The \(2^{k-1}\) products \(Z\) of Section 6 give every solution of \(a^2+b^2=N\) with \(0\le a\le b\) exactly once, with \(a\) the smaller and \(b\) the larger of \(|\mathrm{Re}\,Z|\) and \(|\mathrm{Im}\,Z|\). Consequently \(N\) has exactly \(2^{k-1}\) solutions, and

$$S(N)=\sum_{\eta_2,\dots,\eta_k}\min\bigl(|\mathrm{Re}\,Z|,\,|\mathrm{Im}\,Z|\bigr).$$

Every solution appears. Let \(a^2+b^2=N\) and \(w=a+b\,i\), so that \(w\overline w=N\). The Gaussian prime \(z_{p_1}\) divides \(p_1\), hence \(w\overline w\), hence \(w\) or \(\overline w\); in the second case \(\overline{z_{p_1}}\) divides \(w\). So \(w=\eta_1w_1\) with \(\eta_1\in\{z_{p_1},\overline{z_{p_1}}\}\) and \(|w_1|^2=N/p_1\). Repeating this step with \(p_2,\dots,p_k\) gives \(w=\varepsilon\,\eta_1\eta_2\cdots\eta_k\) with \(|\varepsilon|^2=1\), so \(\varepsilon\) is a unit. If \(\eta_1=\overline{z_{p_1}}\), we replace \(w\) by \(\overline w\), which describes the same solution and has the conjugate factors. Now \(\eta_1=z_{p_1}\), so \(w\) is a unit times one of the products \(Z\), and by Section 4 this \(Z\) describes the solution \((a,b)\).

No solution appears twice. Suppose that two different choices give products \(Z\) and \(Z'\) with the same solution. By Section 4, \(Z'\) is a unit times \(Z\) or a unit times \(\overline Z\). In the first case both sides are products of Gaussian primes with exactly one factor of norm \(p_j\) for every \(j\), so unique factorisation makes each \(\eta_j'\) an associate of \(\eta_j\); since \(z_{p_j}\) and \(\overline{z_{p_j}}\) are not associates, \(\eta_j'=\eta_j\) for every \(j\), so the choices agree. In the second case the factor of norm \(p_1\) is \(z_{p_1}\) on one side and \(\overline{z_{p_1}}\) on the other, which are not associates either. Both cases are impossible.

8. A check with the two-squares counting function

The count agrees with the classical formula \(r_2(n)=4\bigl(d_1(n)-d_3(n)\bigr)\) for the number of ordered pairs \((x,y)\) of integers, signs included, with \(x^2+y^2=n\), where \(d_1(n)\) and \(d_3(n)\) count the divisors of \(n\) that are \(\equiv1\) and \(\equiv3\pmod4\). ^[6] All \(2^k\) divisors of our \(N\) are products of primes \(\equiv1\pmod4\) and therefore \(\equiv1\pmod4\) themselves, so \(r_2(N)=4\cdot2^k\). Each solution accounts for eight of these pairs by Section 4, which leaves \(4\cdot2^k/8=2^{k-1}\) solutions. For \(N=65\) there are \(16\) such pairs.

9. Worked example: the primes 5, 13 and 17

The program arranges the products of Section 6 as a tree. It first fixes the smallest prime of the subset; then, for every larger prime in turn, it leaves the prime out, multiplies by \(z_p\) or multiplies by \(\overline{z_p}\). With \(z_5=1+2i\), \(z_{13}=3+2i\) and \(z_{17}=1+4i\):

• Starting from \(z_5=1+2i\) and leaving out \(13\), the three options for \(17\) give \(1+2i\) for \(N=5\) and \(-7+6i\) and \(9-2i\) for \(N=85\), so \(a=1\), \(6\) and \(2\). Multiplying by \(z_{13}\) gives \(-1+8i\) for \(N=65\), then \(-33+4i\) and \(31+12i\) for \(N=1105\): \(a=1\), \(4\) and \(12\). Multiplying by \(\overline{z_{13}}\) gives \(7+4i\), then \(-9+32i\) and \(23-24i\): \(a=4\), \(9\) and \(23\). These nine leaves add up to \(62\).
• Starting from \(z_{13}=3+2i\), the leaves \(3+2i\), \(-5+14i\) and \(11-10i\) give \(a=2\) for \(N=13\) and \(a=5\) and \(10\) for \(N=221\), together \(17\).
• Starting from \(z_{17}=1+4i\), the single leaf gives \(a=1\).

The \(13\) leaves are exactly the \(13\) solutions of Section 1, and \(62+17+1=80\) is the sum found there. For instance, \((-1+8i)(1+4i)\)\({}=-1-4i+8i-32\)\({}=-33+4i\), and indeed \(33^2+4^2=1105\).

10. Counting the leaves

A number with \(k\) prime factors contributes \(2^{k-1}\) leaves, so by the binomial theorem applied to \((1+2)^{16}\) the whole tree has

$$\begin{aligned}\sum_{k=1}^{16}\binom{16}{k}2^{k-1}&=\frac{3^{16}-1}{2}\\&=21\,523\,360\end{aligned}$$

leaves.

11. How large the numbers get

The largest \(N\), about \(2.5\cdot10^{27}\), is far beyond the range of signed 64-bit integers, whose largest value is about \(9.2\cdot10^{18}\). The program never forms \(N\), however, only the parts of the products. A product \(Z\) of norm \(M\) has \(|\mathrm{Re}\,Z|\le\sqrt M\) and \(|\mathrm{Im}\,Z|\le\sqrt M\), and when \(Z=x+y\,i\) is multiplied by \(z_p\), each of the four partial products such as \(|xu_p|\) is at most \(\sqrt M\sqrt p\). Both bounds are at most the square root of the product of all \(16\) primes, about \(5.0\cdot10^{13}\). For the sum, a solution with \(a\lt b\) has \(2a^2\lt N\), so \(a\lt\sqrt{N/2}\) and \(S(N)\lt2^{k-1}\sqrt{N/2}\), which is \(\prod_{p\mid N}2\sqrt p\) divided by \(2\sqrt2\). Multiplying out \(\prod_p\bigl(1+2\sqrt p\bigr)\) over the \(16\) primes produces one term \(\prod_{p\in T}2\sqrt p\) for every set \(T\) of primes, so adding this bound over all \(N\) gives

$$\begin{aligned}\sum_N S(N)&\lt\frac{\prod_p\bigl(1+2\sqrt p\bigr)-1}{2\sqrt2}\\&\approx3.85\cdot10^{18},\end{aligned}$$

which is below \(2^{63}-1\approx9.22\cdot10^{18}\). So signed 64-bit integers hold every part, every partial product and, as all terms are nonnegative, every partial sum.

How the Code Works

In the C++ version ^[7], main calls solve, which first calls build_representations. Its loops run over odd i and even j with \(i^2+j^2\lt150\) (kLimit) and store Complex(i, j) at reps[i * i + j * j]. For each of the \(16\) primes exactly one such pair exists by Section 5, so reps[p] is \(z_p\); the composite places that the loops also fill, such as \(65\), are never read.

solve then adds up dfs(k + 1, reps[kPrimes[k]]) for \(k=0,1,\dots,15\): kPrimes[k] is the smallest prime of the subset, with its factor fixed to \(z_p\) as in Section 6. In dfs(idx, value), value is the product chosen so far and idx is the index of the next prime to decide. When idx reaches kPrimeCount, which is \(16\), the call returns the smaller absolute value of the two parts, computed with std::llabs and std::min: the \(a\) of Section 7. Otherwise it returns the sum of three calls with idx + 1: with value unchanged, for a prime left out, with value * reps[p] and with value * std::conj(reps[p]). Each leaf is one subset with one choice of factors, hence by Section 7 exactly one solution of one \(N\). At most \(16\) calls of dfs are active at once.

C++ uses std::complex<i64>, where i64 is long long. The C++ standard guarantees std::complex only for floating-point types; for other types the behaviour is unspecified. ^[8] The standard libraries of GCC and Clang both multiply integer complex numbers by the rule of Section 3, which is exact here by Section 11. Java uses a small class Complex with two long fields, whose methods mul and conj return new objects, and the loops of C++ in buildRepresentations. Python builds reps inside solve as tuples of two integers: it tries every odd i and even j below LIMIT and keeps the pairs whose squared sum is below LIMIT, which gives the same prime entries. It multiplies with the helper functions cmul and conj, and its integers have unlimited size. Python and Java return the total from solve as a decimal string, C++ as an integer.

Complexity and Verification

With \(K\) primes the tree has \(L=(3^K-1)/2\) leaves and \(K\) starting calls, and every inner call has three children. Every call except the \(K\) starting ones is the child of exactly one inner call, so the number \(I\) of inner calls satisfies \(3I=L+I-K\), that is, \(I=(L-K)/2\). For \(K=16\) these are \(21\,523\,360\) leaves, one per solution, and \(10\,761\,672\) inner calls, \(32\,285\,032\) calls of dfs in all, and every inner call performs two Gaussian multiplications. Counting each operation on these bounded integers as one step, the time grows like \(3^K\) and the memory is \(O(K)\) for the active calls. Filling the table of \(150\) entries takes \(27\) passes through the inner loop in C++ and Java and \(5550\) in Python. On a desktop computer the C++ and Java programs take well under a second, the Python program about six seconds.

The published programs contain no built-in checks; they print the total. For this article the following was verified separately:

• The tree restricted to the smallest primes gives \(8\) for 5 and 13, \(80\) for 5, 13 and 17, and \(906\), \(12006\) and \(165658\) for the primes below \(30\), \(40\) and \(50\), the same values as a direct search.
• For each of the \(4095\) squarefree numbers built from the first \(12\) primes, 5 to 101, a direct search over all \(a\) with \(2a^2\le N\) finds exactly the solutions produced in Section 6, \(265\,720\) in total.
• A separately written program with 128-bit integers runs through all \(65535\) subsets of the \(16\) primes and their choices of factors with bit masks instead of recursion. For every \(N\) it checks that each product gives \(a^2+b^2=N\) with \(0\lt a\lt b\) and that the \(2^{k-1}\) solutions are pairwise different; the largest part, \(50\,151\,584\,231\,673\), is just below the bound of Section 11. Its total agrees with the output of all three published programs, the C++ one built with GCC and with Clang.

Footnotes and References

The references below supply the background facts used above. The correspondence between choices of conjugate factors and solutions, the size of the tree and the 64-bit bounds are derived explicitly in this article.

1. Project Euler 273 — Sum of Squares. The official statement defines S(N) for the solutions of a² + b² = N with 0 ≤ a ≤ b, gives S(65) = 5 and asks for the sum of S(N) over all squarefree N only divisible by primes of the form 4k + 1 below 150.
2. Wikipedia — Fermat's theorem on sums of two squares. An odd prime p is a sum of two squares if and only if p ≡ 1 (mod 4); such a prime is the norm of a Gaussian prime, p = (x + yi)(x − yi).
3. OEIS A002144 — Pythagorean primes. The primes of the form 4k + 1, starting with the sixteen primes of this problem; a comment notes that their representation as a sum of two nonzero squares is unique.
4. Wikipedia — Gaussian integer. The norm a² + b² is multiplicative, the units are ±1 and ±i, Euclidean division gives unique factorisation up to order and units, and a prime p ≡ 1 (mod 4) is the product of two conjugate, non-associated Gaussian primes.
5. Wikipedia — Brahmagupta–Fibonacci identity. The identity (a² + b²)(c² + d²) = (ac ∓ bd)² + (ad ± bc)² expresses the multiplicativity of the norm a² + b², with its two choices of sign.
6. Wikipedia — Sum of squares function. The function r₂(n) counts the representations n = x² + y² with order and signs, and r₂(n) = 4(d₁(n) − d₃(n)) with the numbers of divisors ≡ 1 and ≡ 3 (mod 4).
7. C++ — Euler273.cpp. The linked immutable C++ revision contains build_representations, dfs and solve. It is the implementation record for this article.
8. cppreference — std::complex. Instantiating std::complex with an integer type for its parts has unspecified behaviour and may fail to compile, so the standard does not guarantee the std::complex<long long> used here.