Detailed mathematical approach

Problem Summary

Let \(F(n)\) be the number of nonzero squarefree Gaussian integers \(z=a+bi\) with norm

$N(z)=a^2+b^2\le n,$

where two values that differ only by multiplication by a unit in \(\{\pm 1,\pm i\}\) are regarded as the same object. The goal is to evaluate \(F(10^{14})\).

A direct scan over all lattice points in the disk \(a^2+b^2\le n\) is far too slow. The solution instead turns squarefreeness into a Möbius-style inclusion-exclusion problem in \(\mathbb{Z}[i]\), collapses that Gaussian sum to an ordinary integer sum indexed by norms, and then accelerates the remaining geometric counts by grouping equal quotients.

Mathematical Approach

Define the lattice-point counting function

$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$

This counts all Gaussian integers of norm at most \(x\), including \(0\). The final formula will use \(Q(x)-1\) so that the origin does not contribute.

Step 1: Write squarefreeness with Gaussian Möbius inversion

In the unique factorization domain \(\mathbb{Z}[i]\), the squarefree indicator has the same shape as over the ordinary integers:

$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d),$

where \(\mu_{\mathbb{Z}[i]}(d)\) is the Gaussian Möbius function. It is \(0\) if \(d\) is divisible by the square of a Gaussian prime, and otherwise \((-1)^k\) when \(d\) is a product of \(k\) distinct Gaussian primes up to units.

Summing over all nonzero Gaussian integers with norm at most \(n\) gives

$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$

The factor \(4\) appears because every nonzero Gaussian integer has exactly four associates obtained by multiplying by the units.

Step 2: Collapse the Gaussian divisor sum by rational norms

Only the norm \(N(d)\) matters in the geometric term, so we group together all Gaussian divisors with the same norm. Define

$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$

Then the previous identity becomes

$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$

This already explains the overall structure of the programs: build the arithmetic weight \(W(m)\), evaluate the circle count \(Q(\cdot)\), and combine them.

Step 3: Derive the prime-power weights from the splitting of rational primes

The function \(W\) is multiplicative, so it is enough to know its values on prime powers. These come directly from the three possible behaviors of a rational prime in \(\mathbb{Z}[i]\).

For \(p=2\), the prime is ramified:

$2=-i(1+i)^2.$

There is only one Gaussian prime above \(2\), with norm \(2\). Therefore the local contribution is

$1-t.$

For \(p\equiv 1\pmod{4}\), the prime splits into two nonassociate conjugate Gaussian primes, both of norm \(p\). The local squarefree possibilities are: choose neither, choose one of the two, or choose both. Hence the local polynomial is

$\left(1-t\right)^2=1-2t+t^2.$

For \(p\equiv 3\pmod{4}\), the prime stays Gaussian-prime, but its norm is \(p^2\). Thus the local polynomial is

$1-t^2.$

Reading off the coefficients gives

$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$

Multiplicativity then determines \(W(m)\) for every \(m\).

Step 4: Count Gaussian integers in a disk

The geometric part is the circle-counting function \(Q(x)\). The implementations use the exact identity

$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$

For each nonnegative \(a\), the inner floor gives the largest admissible \(b>0\). As \(a\) increases, that maximal \(b\) never increases, so a monotone two-pointer scan evaluates the whole sum in \(O(\sqrt{x})\) time rather than recomputing a square root search from scratch for every \(a\).

Step 5: Group equal quotients

The quantity

$k=\left\lfloor\frac{n}{m^2}\right\rfloor$

is constant on an interval of consecutive \(m\)-values. If the current block starts at \(\ell\), its right endpoint is

$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$

If we also precompute the prefix sums

$P(t)=\sum_{m\le t}W(m),$

then the whole block contributes

$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$

This is the key speedup: the expensive circle count is evaluated once per distinct quotient instead of once per \(m\).

Worked Example: \(n=10\)

Here \(\lfloor\sqrt{10}\rfloor=3\), so only \(m=1,2,3\) can appear. From the prime-power rule we get

$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$

Now compute the two needed circle counts:

$Q(10)=37,\qquad Q(2)=9.$

Therefore

$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$

so

$F(10)=7.$

This matches the small checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations follow the same plan. First they set \(M=\lfloor\sqrt{n}\rfloor\) and build a smallest-prime-factor table up to \(M\). This makes it easy to factor every integer \(m\le M\) and to evaluate the multiplicative weight \(W(m)\) from the prime-power rules above.

Next they build prefix sums of \(W\), so any block sum \(\sum_{\ell\le m\le r}W(m)\) is available in constant time. They then partition the interval \(1\le m\le M\) into maximal blocks on which \(\lfloor n/m^2\rfloor\) is constant. For each block they evaluate the lattice count \(Q(k)\), multiply by the corresponding weight sum, add everything together, and divide the final total by \(4\).

The C++ version computes the independent circle counts for different quotient blocks in parallel, while the Python and Java versions perform the same arithmetic sequentially. The mathematical formula is identical in all three languages.

Complexity Analysis

Let \(M=\lfloor\sqrt{n}\rfloor\). Building the smallest-prime-factor data and the multiplicative weights needs \(O(M)\) memory and near-linear time in \(M\). The prefix sums are also \(O(M)\).

The number of distinct values of \(\lfloor n/m^2\rfloor\) is much smaller than \(M\); it is \(O(n^{1/3})\). Each such value requires one circle-count evaluation \(Q(k)\), and that evaluation costs \(O(\sqrt{k})\) time with constant extra memory. In practice the geometric counts dominate the running time, which is why quotient grouping and, in C++, parallel evaluation are the important optimizations.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=556
2. Gaussian integers: Wikipedia — Gaussian integer
3. Gaussian primes and the splitting of rational primes: Wikipedia — Gaussian primes
4. Möbius function: Wikipedia — Möbius function
5. Gauss circle problem: Wikipedia — Gauss circle problem