Solution

Answer: 17971254122360635

Let

$$S(N)=\sum_{n=1}^N s(n),$$

where $s(n)$ is the sum of all Gaussian integer divisors of $n$ whose real part is positive.

We are given:

$$S(10^5)=17924657155.$$

We seek:

$$S(10^8).$$

Mathematical analysis

A Gaussian integer is

$$a+bi,\qquad a,b\in\mathbb Z.$$

Its norm is

$$N(a+bi)=a^2+b^2.$$

A Gaussian integer $a+bi$ divides the rational integer $n$ iff

$$\frac{n}{a+bi}$$

is again a Gaussian integer.

1. Rational divisors ($b=0$)

These are just the ordinary positive divisors of $n$.

Their total contribution is

$$\sum_{n\le N}\sigma(n),$$

where $\sigma(n)$ is the ordinary divisor-sum function.

Using the standard divisor-switching identity:

$$\sum_{n\le N}\sigma(n) = \sum_{d\le N} d\left\lfloor \frac Nd\right\rfloor.$$

2. Non-real Gaussian divisors

Suppose

$$z=a+bi,\qquad b\neq 0.$$

Then

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

For this to be Gaussian integral, both

$$\frac{na}{a^2+b^2}, \qquad \frac{nb}{a^2+b^2}$$

must be integers.

Write

$$a=ga',\qquad b=gb', \qquad \gcd(a',b')=1.$$

Then

$$a^2+b^2=g^2(a'^2+b'^2).$$

Since

$$\gcd(a',a'^2+b'^2)=1, \qquad \gcd(b',a'^2+b'^2)=1,$$

the divisibility condition becomes

$$g(a'^2+b'^2)\mid n.$$

So every Gaussian divisor comes from:

a primitive pair $(a',b')$ with $\gcd(a',b')=1$,
scaled by $g$.

3. Contribution of conjugate pairs

If $a+bi$ divides $n$, then so does $a-bi$.

Their sum is

$$(a+bi)+(a-bi)=2a.$$

So imaginary parts cancel completely.

Thus every primitive pair contributes only real values.

Let

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

For each scale $g$, the divisor pair contributes

$$2ga'$$

to every multiple of $gm$.

Hence total contribution is

$$\sum_{\substack{a',b'>0\ \gcd(a',b')=1}} 2a' \sum_{g\le N/m} g\left\lfloor\frac{N}{gm}\right\rfloor.$$

But

$$\sum_{g\le M} g\left\lfloor\frac{M}{g}\right\rfloor = \sum_{k\le M}\sigma(k).$$

Therefore:

$$S(N) = \sum_{d\le N} d\left\lfloor\frac Nd\right\rfloor + 2\sum_{\substack{a,b>0\ \gcd(a,b)=1}} a, \Sigma!\left(\left\lfloor\frac{N}{a^2+b^2}\right\rfloor\right),$$

where

$$\Sigma(x)=\sum_{n\le x}\sigma(n).$$

This is the key formula.

Efficient algorithm

We need:

Fast computation of

$$\Sigma(x)=\sum_{n\le x}\sigma(n).$$

Using divisor grouping:

$$\Sigma(x) = \sum_{d\le x} d\left\lfloor\frac xd\right\rfloor.$$

This can be evaluated in $O(\sqrt x)$.

Enumerate primitive lattice points

$$a^2+b^2\le N.$$

Only coprime pairs matter.

The resulting algorithm runs comfortably fast in PyPy for $N=10^8$.

Python implementation

from math import gcd, isqrt

N = 10**8

# Summatory sigma:
#   S(n) = sum_{k<=n} sigma(k)
#
# Using:
#   sum sigma(k) = sum d * floor(n/d)
#
# Computed in O(sqrt(n)) using divisor grouping.
def sigma_sum(n):
    total = 0
    d = 1

    while d <= n:
        q = n // d
        r = n // q

        # sum of integers d..r
        block = (d + r) * (r - d + 1) // 2

        total += q * block
        d = r + 1

    return total

# Memoize sigma_sum values because many repeats occur.
memo = {}

def Sigma(n):
    if n not in memo:
        memo[n] = sigma_sum(n)
    return memo[n]

# ----------------------------------------------------
# Rational integer divisors contribution
# ----------------------------------------------------

answer = Sigma(N)

# ----------------------------------------------------
# Non-real Gaussian divisors
# ----------------------------------------------------

limit = isqrt(N)

for a in range(1, limit + 1):

    max_b = isqrt(N - a * a)

    for b in range(1, max_b + 1):

        if gcd(a, b) != 1:
            continue

        m = a * a + b * b

        answer += 2 * a * Sigma(N // m)

print(answer)

Code walkthrough

`sigma_sum(n)`

Computes

$$\sum_{k\le n}\sigma(k).$$

Using

$$\sum_{k\le n}\sigma(k) = \sum_{d\le n} d\left\lfloor\frac nd\right\rfloor.$$

Values of $\lfloor n/d\rfloor$ stay constant on intervals, so we group them in $O(\sqrt n)$.

Rational divisor contribution

answer = Sigma(N)

This accounts for all ordinary positive divisors.

Primitive Gaussian divisors

We loop over coprime pairs:

if gcd(a, b) != 1:
    continue

with

m = a*a + b*b

Every such primitive Gaussian divisor contributes

$$2a\cdot \Sigma!\left(\left\lfloor\frac{N}{m}\right\rfloor\right).$$

So:

answer += 2 * a * Sigma(N // m)

Verification on small cases

For $N=5$:

$$S(5)=35,$$

matching the statement.

For $N=10^5$:

$$S(10^5)=17924657155,$$

also matching the statement.

So the derivation and implementation are correct.

Final computation

Evaluating the program for

$$N=10^8$$

gives:

$$17971254122360635$$

Answer: 17971254122360635

Project Euler Problem 153