Detailed mathematical approach

Problem Summary

Let \(s(n)\) denote the Gaussian-divisor sum defined in the problem, and let

$S(N)=\sum_{n=1}^{N} s(n),\qquad N=10^8.$

Ordinary positive divisors of \(n\) are also Gaussian divisors, purely imaginary divisors have real part \(0\), and genuinely complex divisors occur in conjugate pairs whose imaginary parts cancel. The task is therefore not to factor every \(n\) separately in \(\mathbb{Z}[i]\), but to count how often each possible Gaussian divisor appears among the integers up to \(N\).

The C++, Python, and Java implementations all exploit the same reorganization: isolate the ordinary divisor contribution, classify every non-real Gaussian divisor by a primitive direction and an integer scale, and evaluate the resulting sums through a fast summatory divisor routine.

Mathematical Approach

The entire solution comes from turning a sum over integers \(n\) into a sum over Gaussian divisors that can divide some \(n\le N\).

Ordinary divisors versus genuine Gaussian divisors

The divisors on the real axis are exactly the ordinary positive divisors. Their total contribution to \(S(N)\) is

$T(N)=\sum_{m=1}^{N}\sigma(m)=\sum_{d=1}^{N} d\left\lfloor\frac{N}{d}\right\rfloor,$

because each positive integer \(d\) divides exactly \(\lfloor N/d\rfloor\) integers up to \(N\). This is the familiar summatory divisor function.

Purely imaginary divisors need no separate treatment: their real part is zero, so they do not affect the final total. The only additional work comes from Gaussian divisors with both real and imaginary parts nonzero.

Primitive directions and scaled divisors

Every contributing non-real Gaussian divisor can be represented uniquely in the form

\(d(a+bi)\) or \(d(a-bi)\).

where \(d\ge 1\), \(a,b\ge 1\), and \(\gcd(a,b)=1\). Any common factor of the real and imaginary parts is absorbed into the positive integer scale \(d\), so \((a,b)\) describes a primitive Gaussian direction.

For that primitive direction, define the norm

$m=a^2+b^2.$

This is exactly why the implementations loop only over coprime pairs \((a,b)\): non-primitive pairs would merely duplicate divisors already counted by a larger scale \(d\).

The divisibility test for ordinary integers

Fix a primitive direction \((a,b)\) and a scale \(d\). To determine when \(d(a+bi)\) divides an ordinary integer \(n\), multiply by the conjugate:

$\frac{n}{d(a+bi)}=\frac{n(a-bi)}{d(a^2+b^2)}=\frac{n(a-bi)}{dm}.$

Because \(\gcd(a,b)=1\), we also have \(\gcd(a,m)=\gcd(b,m)=1\). Therefore the real and imaginary parts of the quotient are integers exactly when \(dm\mid n\).

This is the crucial simplification. Divisibility by a Gaussian number is reduced to an ordinary divisibility condition involving the scale \(d\) and the primitive norm \(m=a^2+b^2\).

Collapsing the global sum

For fixed \((a,b)\) and \(d\), the conjugate pair \(d(a+bi)\) and \(d(a-bi)\) contributes the same real part \(da\), so together they contribute \(2da\) to every integer \(n\) divisible by \(dm\). There are exactly \(\left\lfloor N/(dm)\right\rfloor\) such integers up to \(N\).

Hence one primitive direction \((a,b)\) contributes

$\sum_{d\ge 1} 2da\left\lfloor\frac{N}{dm}\right\rfloor=2a\sum_{d\le N/m} d\left\lfloor\frac{N/m}{d}\right\rfloor=2a\,T\!\left(\left\lfloor\frac{N}{m}\right\rfloor\right).$

Adding the ordinary divisor part gives the exact formula implemented by all three languages:

$\boxed{S(N)=T(N)+\sum_{\substack{a,b\ge 1\\ \gcd(a,b)=1\\ a^2+b^2\le N}} 2a\,T\!\left(\left\lfloor\frac{N}{a^2+b^2}\right\rfloor\right).}$

The condition \(a^2+b^2\le N\) is automatic: if the primitive norm already exceeds \(N\), then even the unscaled Gaussian divisor cannot divide any integer at most \(N\).

Worked example: the primitive pair \(1+2i\)

Take \(a=1\), \(b=2\), so \(m=1^2+2^2=5\). If \(N=20\), then the relevant scaled divisors are

$1\pm 2i,\quad 2\pm 4i,\quad 3\pm 6i,\quad 4\pm 8i,$

because \(d\le N/m=4\). The pair \(d(1\pm 2i)\) divides exactly the multiples of \(5d\), and its total real contribution is \(2d\cdot \lfloor 20/(5d)\rfloor\). Summing over \(d=1,2,3,4\) gives

$2\bigl(1\cdot 4+2\cdot 2+3\cdot 1+4\cdot 1\bigr)=30.$

This is precisely

$2a\,T\!\left(\left\lfloor\frac{20}{5}\right\rfloor\right)=2\,T(4),$

which shows on a small example how one primitive direction collapses to one summatory-divisor query.

How the Code Works

Fast evaluation of the divisor-prefix function

The first ingredient is \(T(M)=\sum_{d\le M} d\lfloor M/d\rfloor\). The implementation does not iterate over all \(d\) one by one. Instead it groups together intervals on which the quotient \(\lfloor M/d\rfloor\) is constant.

If the current left endpoint is \(l\), then \(q=\lfloor M/l\rfloor\) stays unchanged up to \(r=\lfloor M/q\rfloor\). The whole block contributes

$q\sum_{d=l}^{r} d=q\cdot \frac{(l+r)(r-l+1)}{2}.$

This standard equal-quotient grouping reduces one prefix query from linear time to about \(O(\sqrt M)\).

Enumerating primitive Gaussian directions

The implementations start the answer with the ordinary divisor contribution \(T(N)\). They then scan all integer pairs \(a\ge 1\), \(b\ge 1\) satisfying \(a^2+b^2\le N\), and retain only the coprime pairs. Those are the primitive directions from the derivation above.

For each surviving pair, the norm \(a^2+b^2\) is computed, the value \(T\!\left(\left\lfloor N/(a^2+b^2)\right\rfloor\right)\) is queried, multiplied by \(2a\), and added to the running total.

Avoiding repeated work

Many different primitive pairs produce the same quotient \(\left\lfloor N/(a^2+b^2)\right\rfloor\). The C++ and Python implementations store previously computed divisor-prefix values so that repeated queries are answered immediately. The Java implementation applies the same formula directly without that memoization layer, but the mathematical accumulation is identical in all three languages.

Complexity Analysis

The dominant scan is the enumeration of lattice points in the quarter-disk \(a^2+b^2\le N\), which contains \(O(N)\) integer pairs because the radius is \(\sqrt N\). The coprimality test changes only the constant factor.

Each distinct divisor-prefix query is evaluated in about \(O(\sqrt M)\) time by equal-quotient grouping. In the memoized variants, many queries reuse an earlier result because the same quotient occurs repeatedly. Memory usage is \(O(U)\) for the number \(U\) of cached prefix values; aside from that cache, the algorithm uses only constant extra loop state.

The real gain is structural rather than cosmetic: instead of factoring every integer up to \(10^8\) in \(\mathbb{Z}[i]\), the solution performs one global scan of primitive Gaussian directions and a manageable number of summatory-divisor evaluations.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=153
2. Gaussian integers: Wikipedia - Gaussian integer
3. Coprime integers: Wikipedia - Coprime integers
4. Divisor function: Wikipedia - Divisor function
5. Divisor summatory function: Wikipedia - Divisor summatory function