Detailed mathematical approach

Problem Summary

Let

$P=\{5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149\}$

be the primes below \(150\) with \(p\equiv1\pmod4\). For every non-empty squarefree product of distinct primes from \(P\), the problem considers all representations

$n=a^2+b^2$

and asks us to add the smaller component from each representation. The final numeric answer is omitted here; the goal is to explain why the code can enumerate all valid representations without overcounting.

Mathematical Approach

The whole solution is built on the arithmetic of Gaussian integers \(\mathbb Z[i]\), where the norm

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

is multiplicative.

1) Each eligible prime splits as a sum of two squares

By Fermat's two-square theorem, every prime \(p\equiv1\pmod4\) has a representation

$p=u_p^2+v_p^2.$

Equivalently, in \(\mathbb Z[i]\) it factors as

$p=(u_p+iv_p)(u_p-iv_p)=z_p\overline{z_p}.$

The code stores one canonical factor \(z_p\) for each prime. Its helper build_representations() searches odd \(u\) and even \(v\), which is enough for these primes. For example,

$5=1^2+2^2,\qquad 13=3^2+2^2,\qquad 17=1^2+4^2.$

So the stored Gaussian representatives are \(z_5=1+2i\), \(z_{13}=3+2i\), \(z_{17}=1+4i\).

2) Multiplication gives the sum-of-two-squares composition law

If \(z=a+bi\) and \(w=c+di\), then

$zw=(ac-bd)+i(ad+bc),$

and therefore

$N(zw)=N(z)N(w).$

Written out in ordinary integers, this is the classical identity

$\left(a^2+b^2\right)\left(c^2+d^2\right)=(ac-bd)^2+(ad+bc)^2.$

If we conjugate one factor we get the second standard form

$\left(a^2+b^2\right)\left(c^2+d^2\right)=(ac+bd)^2+(ad-bc)^2.$

This is exactly why multiplying Gaussian factors generates every representation of a product as a sum of two squares.

3) All representations of a squarefree product come from conjugate choices

Take a non-empty subset \(S\subseteq P\), and define

$n=\prod_{p\in S} p.$

Since every such prime splits as \(p=z_p\overline{z_p}\), a factorization of \(n\) in \(\mathbb Z[i]\) is obtained by choosing exactly one factor from each pair:

$Z=\prod_{p\in S}\eta_p,\qquad \eta_p\in\{z_p,\overline{z_p}\}.$

Then \(N(Z)=n\). If \(Z=x+iy\), we immediately get

$n=x^2+y^2.$

Conversely, for a squarefree product of such primes, every representation arises from such a choice, up to the usual symmetries \((x,y)\mapsto(\pm x,\pm y)\) and swapping the two coordinates.

4) Why the DFS has exactly three branches per later prime

Suppose we have already fixed some current Gaussian product \(Z\). For the next prime \(p\), there are only three meaningful options:

1. skip \(p\), meaning \(p\) is not in the subset;

2. multiply by \(z_p\);

3. multiply by \(\overline{z_p}\).

That is why the recursion is ternary:

$\text{dfs}(idx,Z)=\text{dfs}(idx+1,Z)+\text{dfs}(idx+1,Zz_p)+\text{dfs}(idx+1,Z\overline{z_p}).$

At a leaf, the code has built one complete representation \(Z=x+iy\), so it contributes

$\min(|x|,|y|).$

5) Why the outer loop fixes the first chosen prime

A naive DFS starting from \(1\) and allowing every selected prime to choose either \(z_p\) or \(\overline{z_p}\) would double-count every non-empty subset. The reason is global conjugation:

$Z=\prod_{p\in S}\eta_p \qquad \Longrightarrow \qquad \overline{Z}=\prod_{p\in S}\overline{\eta_p}.$

But \(Z=x+iy\) and \(\overline{Z}=x-iy\) represent the same unordered pair \((|x|,|y|)\), hence the same contribution \(\min(|x|,|y|)\).

The implementation removes that symmetry by choosing the first included prime \(p_k\) in the outer loop and forcing its factor to be \(z_{p_k}\), not \(\overline{z_{p_k}}\). Only primes with larger index branch between \(z_p\) and \(\overline{z_p}\). So each non-empty subset is enumerated exactly once up to global conjugation.

6) Worked example: \(65=5\cdot13\)

Use \(z_5=1+2i\) and \(z_{13}=3+2i\). The two conjugate choices for the second prime produce the two relevant representations:

$\left(1+2i\right)\left(3+2i\right)=-1+8i \quad\Longrightarrow\quad 65=1^2+8^2,$

$\left(1+2i\right)\left(3-2i\right)=7+4i \quad\Longrightarrow\quad 65=7^2+4^2.$

Therefore the total contribution of \(65\) is

$\min(1,8)+\min(7,4)=1+4=5.$

This also gives a clean checkpoint for the first two primes only:

$5\mapsto1,\qquad 13\mapsto2,\qquad 65\mapsto5,\qquad \text{total}=8.$

The same code, restricted to the first three primes \(\{5,13,17\}\), gives \(80\), and restricted to the first four gives \(906\).

How the Code Works

Representation table. build_representations() scans odd \(i\) and even \(j\), and when \(i^2+j^2<150\) it stores the Gaussian integer \(i+ji\). For the 16 target primes, this yields exactly one canonical representation per prime.

Recursive enumeration. dfs(idx, value) processes the remaining primes in increasing order. The three recursive calls correspond to skip, multiply by \(z_p\), and multiply by \(\overline{z_p}\).

Leaf contribution. When all primes have been processed, the current Gaussian product value = x + iy encodes one valid representation of one squarefree product. The code returns min(abs(x), abs(y)).

Canonical start. solve() loops over the position of the first chosen prime and starts the DFS with reps[p]. This is the exact mechanism that removes the global conjugation duplicate.

Complexity Analysis

If \(k=16\) is the number of eligible primes, the direct search has size

$\sum_{m=0}^{k-1}3^m=\frac{3^k-1}{2}=O(3^k).$

That sounds large in theory, but here \(k\) is tiny and each state only performs a small Gaussian-integer multiplication on 64-bit integers. The recursion depth is \(O(k)\), so memory usage is \(O(k)\).

Further Reading

1. Problem page: https://projecteuler.net/problem=273
2. Fermat's theorem on sums of two squares: https://en.wikipedia.org/wiki/Fermat%27s_theorem_on_sums_of_two_squares
3. Gaussian integers and norms: https://en.wikipedia.org/wiki/Gaussian_integer
4. Brahmagupta-Fibonacci two-square identity: https://en.wikipedia.org/wiki/Brahmagupta%E2%80%93Fibonacci_identity