Detailed mathematical approach

Problem Summary

We must count integers \(n \le 2\times 10^9\) that can all be written in the four forms

$n=a^2+b^2=a^2+2b^2=a^2+3b^2=a^2+7b^2,$

with positive integers \(a\) and \(b\) in each representation.

The implementations do not search directly over all pairs \((a,b)\). Instead they first count a slightly larger set in which zero is temporarily allowed, classify those integers by their squarefree kernel, and then repair the overcount caused by perfect squares.

Mathematical Approach

Let \(L=2\times 10^9\) and let \(D=\{1,2,3,7\}\). For the main counting argument it is convenient to relax the condition and first consider integers \(n\) such that for every \(d\in D\) there exist \(a_d,b_d\ge 0\) with

$n=a_d^2+d\,b_d^2.$

The strict positivity condition will be restored afterward.

Squarefree Kernels Control the Relaxed Problem

Every positive integer has a unique decomposition

$n=q\,m^2,$

where \(q\) is squarefree. The implementations take as their number-theoretic starting point the characterization that simultaneous representability in all four relaxed forms depends only on this squarefree kernel. More precisely, the relaxed condition holds exactly when every prime factor of \(q\) lies in the residue classes

$p\equiv 1,\ 25,\ 121 \pmod{168}.$

So the admissible kernels are precisely

$q=\prod_{i=1}^{r} p_i,$

where the \(p_i\) are distinct primes from those three classes; \(q=1\) is the empty product. Once such a kernel is fixed, every number of the form \(q\,m^2\) belongs to the relaxed set, and once the kernel contains any other prime, the number can be discarded immediately.

Counting Admissible Kernels

For a fixed admissible kernel \(q\), the valid integers are exactly the numbers \(q\,m^2\le L\). Their count is therefore

$\left\lfloor \sqrt{\frac{L}{q}} \right\rfloor.$

If \(R(L)\) denotes the relaxed count, then

$R(L)=\sum_{\substack{q\le L\\ q\text{ admissible}}}\left\lfloor \sqrt{\frac{L}{q}} \right\rfloor.$

The implementations split this sum into three disjoint cases.

The kernel \(q=1\) contributes \(\lfloor\sqrt{L}\rfloor\). A singleton kernel \(q=p\) contributes \(\lfloor\sqrt{L/p}\rfloor\) for each good prime \(p\). All remaining kernels contain at least two distinct good primes, so the implementations enumerate their products by depth-first search and add the same floor term for each feasible product.

An important pruning rule comes from the smallest good prime. Since the first prime with residue \(1\), \(25\), or \(121\) modulo \(168\) is \(193\), any admissible kernel with at least two prime factors is at least \(193^2\). That is why the multi-prime stage can stop very early for small partial products.

Why Perfect Squares Need a Separate Correction

The relaxed count is not yet the answer, because the original problem requires positive \(a\) and \(b\), not merely nonnegative ones. The only systematic overcount comes from perfect squares.

Indeed, every square \(m^2\) is automatically included in the relaxed set via the trivial identities

$m^2=m^2+d\cdot 0^2 \qquad (d\in D).$

For a nonsquare, the degenerate case \(b=0\) cannot occur. And if \(a=0\) for \(d=2\), \(3\), or \(7\), then the squarefree kernel would contain \(2\), \(3\), or \(7\), which is impossible for an admissible kernel. So the positivity repair can be written as

$A(L)=R(L)-\lfloor\sqrt{L}\rfloor+G(L),$

where \(G(L)\) counts the squares \(m^2\le L\) that really do admit positive representations for all four coefficients.

Turning Square Representations into a Divisor Problem

To compute \(G(L)\), fix \(d\in D\) and ask when a square satisfies

$m^2=a^2+d\,b^2$

with \(a,b>0\). Rearranging gives

$(m-a)(m+a)=d\,b^2.$

If we set

$u=m-a,\qquad v=m+a,$

then we obtain

$uv=d\,b^2,\qquad u<v,\qquad u\equiv v\pmod 2,$

and conversely any factor pair \((u,v)\) satisfying these conditions yields

$m=\frac{u+v}{2},\qquad a=\frac{v-u}{2}.$

That turns the search for positive square representations into a divisor enumeration problem. For each fixed \(d\) and each positive \(b\), factor \(d\,b^2\), generate all divisors \(u\), define \(v=(d\,b^2)/u\), keep only pairs with the same parity and \(v>u\), and mark the resulting \(m\). Doing this separately for \(d=1,2,3,7\) produces four boolean tables over \(m\le \sqrt{L}\), and their intersection is exactly the set counted by \(G(L)\).

Worked Example: The Square \(3600\)

The number \(3600=60^2\) illustrates why the correction term is necessary. Its squarefree kernel is \(q=1\), so it is counted automatically in the relaxed total \(R(L)\).

It also survives the positivity correction, because it genuinely has positive representations in all four forms:

$3600=48^2+36^2=20^2+2\cdot 40^2=30^2+3\cdot 30^2=45^2+7\cdot 15^2.$

So \(3600\) is first counted among all squares through the kernel \(q=1\), then removed by the blanket subtraction of \(\lfloor\sqrt{L}\rfloor\), and finally added back because it belongs to \(G(L)\).

How the Code Works

The C++, Python, and Java implementations follow exactly the split described above. First they compute the relaxed quantity \(R(L)\). The singleton-kernel contribution is obtained by scanning primes up to \(L\) with a segmented sieve and keeping only primes in the three admissible residue classes modulo \(168\). Each such prime \(p\) contributes \(\left\lfloor \sqrt{L/p} \right\rfloor\).

The multi-prime part uses the ordered list of good primes up to \(L/193\). A depth-first search grows products of distinct good primes in increasing order, stops as soon as the next multiplication would exceed \(L\), and adds \(\left\lfloor \sqrt{L/q} \right\rfloor\) whenever the current product \(q\) already contains at least two primes. In this way every admissible multi-prime kernel is visited exactly once.

The square correction then builds a smallest-prime-factor table up to \(\sqrt{L}\). For each \(d\in\{1,2,3,7\}\) and each positive \(b\le \sqrt{L/d}\), the implementation factors \(d\,b^2\), generates all divisors, reconstructs every valid \(m\) from the factor pairs \((u,v)\), and marks it. The final answer is the relaxed kernel count minus all squares plus the number of \(m\) marked in every one of the four tables.

The language-specific execution strategy differs slightly, but the mathematics is the same. The C++ and Java implementations parallelize the segmented prime sweep across workers, while the Python implementation uses multiple processes for that same stage.

Complexity Analysis

The dominant large-scale cost is the singleton-kernel stage, which performs a segmented prime sweep over the interval up to \(L\) while sieving with the base primes up to \(\sqrt{L}\). In asymptotic terms this is standard segmented-sieve work over a range of length \(L\), while memory stays bounded by the segment size and the small base-prime list.

The multi-prime DFS is much smaller in practice, because admissible kernels grow quickly once distinct good primes are multiplied together, and the recursion prunes immediately when the running product would exceed \(L\). The square-correction phase only works up to \(\sqrt{L}\), using a smallest-prime-factor table and divisor generation for numbers of the form \(d\,b^2\).

Memory usage is modest throughout: one segment buffer for the sieve, prime-factor data up to \(\sqrt{L}\), and four mark arrays of length \(\lfloor\sqrt{L}\rfloor\). The whole algorithm is efficient because it replaces a huge search in the \((a,b)\)-plane by kernel counting plus a focused correction on squares.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=229
2. Binary quadratic form: Wikipedia - Binary quadratic form
3. Sum of two squares theorem: Wikipedia - Sum of two squares theorem
4. Squarefree integer: Wikipedia - Squarefree integer
5. Segmented sieve: cp-algorithms - Sieve of Eratosthenes