Problem 153: Investigating Gaussian Integers

Project Euler Problem 153 Solution

EulerSolve provides an optimized solution for Project Euler Problem 153, Investigating Gaussian Integers, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

Detailed mathematical approach

Problem Summary

A Gaussian integer is a complex number $x+yi$, where $i^2=-1$, whose real part $x$ and imaginary part $y$ are both integers. It is called a divisor of a rational (ordinary) integer $n$ if the quotient $n/(x+yi)$ is again a Gaussian integer. For $n\ge1$, let $s(n)$ be the sum of the divisors of $n$ whose real part $x$ is positive; divisors such as $i$ or $2i$, with real part $0$, are not counted. For example, $5$ has the six divisors $1$, $1+2i$, $1-2i$, $2+i$, $2-i$ and $5$ with positive real part, so $s(5)=12$, and the table of the statement gives $s(1)=1$, $s(2)=5$, $s(3)=4$, $s(4)=13$ and $s(5)=12$. With

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

the statement gives $S(5)=35$ and $S(10^5)=17924657155$, and we must find $S(10^8)$. ^[1]

The whole solution in one sentence: every divisor with positive real part is an ordinary divisor or a multiple $d(a\pm bi)$ of a Gaussian integer with coprime $a,b\ge1$, such a multiple divides $n$ exactly when the ordinary number $d(a^2+b^2)$ does, and so the total turns into ordinary divisor sums with simple weights, which are added up quickly in blocks of equal quotients.

Mathematical Approach

Instead of finding the divisors of each $n\le10^8$ separately, we describe all of them by an ordinary divisibility condition and then exchange the order of summation, so that each possible divisor is handled once for all $n\le N$ together.

1. Testing a divisor by hand

The statement divides $5$ by $1+2i$ by multiplying numerator and denominator by the complex conjugate $1-2i$. In general the conjugate of $z=x+yi$ is $\bar z=x-yi$, and $z\bar z=x^2+y^2$ is the norm of $z$. ^[2] For $z\ne0$ the same trick gives

$$\begin{aligned}\frac{n}{x+yi}&=\frac{n(x-yi)}{x^2+y^2}\\&=\frac{nx}{x^2+y^2}-\frac{ny}{x^2+y^2}\,i,\end{aligned}$$

so $x+yi$ divides $n$ exactly when its norm divides both $nx$ and $ny$. For $n=2$ and $z=1+i$, the norm $2$ divides both $nx=2$ and $ny=2$, and the quotient is $1-i$. A divisor is never larger than $n$ in absolute value, because $n=zw$ with a nonzero Gaussian integer $w$, and $|w|\ge1$. Testing the finitely many candidates with $x\gt0$ and $x^2+y^2\le n^2$ gives the four divisors $1$, $1+i$, $1-i$ and $2$ of $2$, so $s(2)=5$.

Two facts follow. First, if $z$ divides $n$, so does $\bar z$, because the conjugate of the quotient $n/z$ is $n/\bar z$. The divisors with $y\ne0$ therefore come in pairs $x\pm yi$ whose imaginary parts cancel, and $s(n)$ is the sum of the real parts, an ordinary positive integer. Second, a real divisor $x\gt0$ leaves the real quotient $n/x$, which is a Gaussian integer only when it is an ordinary integer. So the real divisors are the ordinary positive divisors of $n$, and their sum is $\sigma(n)$, the classical sum of divisors. For $n=5$ these are $1$ and $5$, the pairs $1\pm2i$ and $2\pm i$ add $2\cdot1+2\cdot2=6$, and $s(5)=6+6=12$.

2. Primitive pairs

Take a divisor $x+yi$ with $x\gt0$ and $y\ne0$ and pull out the greatest common divisor of its parts: with $d=\gcd(x,|y|)$, $a=x/d$ and $b=|y|/d$, it equals $d(a+bi)$ if $y\gt0$ and $d(a-bi)$ if $y\lt0$, where $a,b\ge1$ and $\gcd(a,b)=1$. This form is unique, because $\gcd(x,|y|)=d\gcd(a,b)$ forces $d$. For example, $2+2i=2(1+i)$ and $3-6i=3(1-2i)$. Geometrically, $a\pm bi$ is the first lattice point on the ray from $0$ through $x+yi$. We call $(a,b)$ a primitive pair and $m=a^2+b^2$ its norm; the divisor $d(a\pm bi)$ has the real part $da$.

3. When a multiple of a primitive pair divides $n$

Claim. Let $a,b\ge1$ be coprime, $m=a^2+b^2$, and $d,n\ge1$. Then $d(a+bi)$ divides $n$ if and only if $dm$ divides $n$, and the same holds for $d(a-bi)$.

Proof. The norm of $d(a+bi)$ is $d^2m$, so by Section 1 the quotient is

$$\frac{n}{d(a+bi)}=\frac{na}{dm}-\frac{nb}{dm}\,i.$$

If $dm$ divides $n$, both parts are integers. Conversely, let $dm$ divide $na$ and $nb$. Since $\gcd(a,b)=1$, Bézout's identity gives integers $u$ and $v$ with $ua+vb=1$, ^[3] so $dm$ also divides $u\cdot na+v\cdot nb=n$. Finally, $n/(d(a-bi))$ is the conjugate of the quotient above, so it is a Gaussian integer exactly when that one is. This proves the claim.

The claim reproduces the statement: $1\pm2i$ and $2\pm i$ have $m=5$ and divide $5$, whereas $1+i$ has $m=2$, and $2$ does not divide $5$. Coprimality matters: $2+2i$ divides $4$, with quotient $1-i$, although its norm $8$ does not divide $4$; written as $2(1+i)$ it has $dm=2\cdot2=4$, which does. The condition depends on $a$ and $b$ only through $m$, so the four primitive pairs $(1,8)$, $(4,7)$, $(7,4)$ and $(8,1)$ of norm $65$ give divisors of exactly the same numbers. A divisor $d(a\pm bi)$ of some $n\le N$ needs $dm\le N$, so only pairs with $m\le N$ matter.

4. Summing over the divisors instead of over $n$

By Sections 1 to 3, with $(a,b)$ running over the primitive pairs and the factor $2$ counting both conjugates,

$$s(n)=\sigma(n)+\sum_{(a,b)}\ \sum_{d\ge1,\ dm\mid n}2da.$$

We add this over $n=1,\dots,N$ and exchange the order of summation, using that a positive integer $k$ divides exactly $\lfloor N/k\rfloor$ of the numbers $1,\dots,N$. The ordinary divisors give

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

the summatory sum-of-divisors function, with $T(1)=1$, $T(2)=4$, $T(5)=21$ and $T(10)=87$. ^[4] For a primitive pair $(a,b)$ and a multiplier $d$, the two divisors $d(a\pm bi)$ add $2da$ to $s(n)$ for each of the $\lfloor N/(dm)\rfloor$ multiples $n$ of $dm$ up to $N$. These multiples are the numbers $mk$ with $k\le\lfloor N/m\rfloor$ and $d\mid k$, so with $M=\lfloor N/m\rfloor$ their number is $\lfloor M/d\rfloor$, an instance of the rule for nested divisions. ^[5] The pair therefore contributes $2a\,T(M)$, since $\sum_d d\lfloor M/d\rfloor=T(M)$, and

$$\begin{aligned}S(N)&=T(N)\\&+\sum_{\substack{a,b\ge1\\ \gcd(a,b)=1\\ a^2+b^2\le N}}2a\,T\Bigl(\Bigl\lfloor\frac{N}{a^2+b^2}\Bigr\rfloor\Bigr).\end{aligned}$$

5. Worked example: the table of the statement

For $N=5$ we need $T(5)=1+3+4+7+6=21$ and the primitive pairs with $a^2+b^2\le5$, namely $(1,1)$ with $m=2$, and $(1,2)$ and $(2,1)$ with $m=5$:

$(1,1)$ gives $2\cdot1\cdot T(2)=8$: in the table, $1\pm i$ divide $2$ and $4$, with real parts $1+1$ each time, and $2\pm2i$ divide $4$, with real parts $2+2$;
$(1,2)$ gives $2\cdot1\cdot T(1)=2$, from $1\pm2i$ in $n=5$;
$(2,1)$ gives $2\cdot2\cdot T(1)=4$, from $2\pm i$ in $n=5$.

Hence $S(5)=21+8+2+4=35$, as stated. For $N=10$ the pairs $(1,3)$ and $(3,1)$ of norm $10$ join, and $\lfloor10/2\rfloor=5$ and $\lfloor10/5\rfloor=2$. Grouping the pairs of equal norm,

$$\begin{aligned}S(10)&=T(10)+2T(5)\\&+(2+4)\,T(2)+(2+6)\,T(1)\\&=87+42+24+8=161.\end{aligned}$$

The pair $(2,2)$ is not primitive and does not occur: the divisors $2\pm2i$ are already counted with $d=2$ in the pair $(1,1)$.

6. Evaluating $T(M)$ in blocks

The sum $T(M)=\sum_{l\le M}l\lfloor M/l\rfloor$ has $M$ terms, but its quotients change rarely: for $M=10$ the values $\lfloor10/l\rfloor$ are $10,5,3,2,2,1,1,1,1,1$. Precisely, $\lfloor M/j\rfloor=q$ holds exactly when $q\le M/j\lt q+1$, that is, when $M/(q+1)\lt j\le M/q$. So the block of equal quotients that starts at $l$, with $q=\lfloor M/l\rfloor$, ends at $r=\lfloor M/q\rfloor$ and contributes

$$\begin{aligned}&q\,\bigl(l+(l+1)+\dots+r\bigr)\\&\quad=q\cdot\frac{(l+r)(r-l+1)}{2}.\end{aligned}$$

There are at most $2\sqrt M$ blocks: at most $\sqrt M$ of them start at some $l\le\sqrt M$, and a block that starts at $l\gt\sqrt M$ has a quotient $q\le M/l\lt\sqrt M$, different for different blocks. For $M=10$ the blocks are $\{1\}$, $\{2\}$, $\{3\}$, $\{4,5\}$ and $\{6,\dots,10\}$, and

$$\begin{aligned}T(10)&=10\cdot1+5\cdot2+3\cdot3\\&\quad+2\cdot9+1\cdot40=87.\end{aligned}$$

For $M=10^8$ there are $19\,999$ blocks instead of $10^8$ terms. A related idea, Dirichlet's hyperbola method, evaluates the divisor summatory function $\sum_{k\le x}\lfloor x/k\rfloor$ in $O(\sqrt x)$ time by splitting the lattice points under the hyperbola $jk=x$ at $\sqrt x$. ^[6]

7. Few different arguments

The formula of Section 4 asks for $T(\lfloor N/m\rfloor)$ once for every primitive pair, but by the argument of Section 6 the numbers $\lfloor N/m\rfloor$ take at most $2\sqrt N$ different values. Storing every computed $T(M)$ in a table and looking it up when the same $M$ comes again therefore computes each value only once. For $N=10^8$ only $7\,438$ different arguments occur among the $47\,746\,495$ primitive pairs and $N$ itself, fewer than $2\sqrt N=20\,000$, because $m$ runs only over the norms of primitive pairs.

8. How large the total can be

Every term of $T(M)$ satisfies $l\lfloor M/l\rfloor\le M$, so $T(M)\le M^2$, and a primitive pair contributes at most $2aN^2/(a^2+b^2)^2$. For fixed $a$ put $f(t)=1/(a^2+t^2)^2$. Since $f$ decreases, each term $f(b)$ is at most the area under $f$ between $b-1$ and $b$, so the sum over $b$ is at most the integral below. The substitution $b=a\tan\theta$ turns that integral into the integral of $\cos^2\theta/a^3$ over $0\le\theta\le\pi/2$, and since $\cos^2\theta$ and $\sin^2\theta$ have the same integral over this range and add up to $1$, that integral is $\pi/4$ times $1/a^3$:

$$\begin{aligned}\sum_{b\ge1}\frac{1}{(a^2+b^2)^2}&\le\int_0^\infty\frac{db}{(a^2+b^2)^2}\\&=\frac{\pi}{4a^3}.\end{aligned}$$

Together with the telescoping bound

$$\sum_{a\ge1}\frac{1}{a^2}\le1+\sum_{a\ge2}\Bigl(\frac{1}{a-1}-\frac{1}{a}\Bigr)=2$$

this gives

$$\begin{aligned}S(N)&\le N^2+\sum_{a\ge1}2a\cdot\frac{\pi N^2}{4a^3}\\&\le(1+\pi)\,N^2,\end{aligned}$$

which is below $4.2\cdot10^{16}$ for $N=10^8$, far below $2^{63}\approx9.2\cdot10^{18}$. Signed and unsigned 64-bit integers therefore hold every partial sum.

How the Code Works

The C++ program ^[7] keeps all sums in u64, an alias of std::uint64_t. Its main first reads the options with parse_arguments: --limit=N replaces the default $N=10^8$ and is read by parse_int_after_prefix, which accepts only a nonempty string of decimal digits but accumulates it in an int without a range check, so an out-of-range value causes signed overflow, undefined behaviour in C++; and --skip-checkpoints turns the checks off. An unknown argument or a limit below $1$ ends the program with exit status $1$. Otherwise run_checkpoints compares solve(5) with $35$ and solve(100000) with $17924657155$, the two values of the statement, and the program exits with status $2$ if either differs. Finally main prints solve(limit).

solve implements the formula of Section 4. It creates the table of Section 7, a std::unordered_map<int, u64> with room reserved for $32\,768$ entries, and starts with answer = sigma_prefix(limit, …), which is $T(N)$. The outer loop runs $a=1,2,\dots$ while $a^2\le N$, and bmax, the largest $b$ with $a^2+b^2\le N$, is the square root of $N-a^2$ in double precision, truncated to an integer; for $N=10^8$ a check over all $a$ confirmed that this equals $\lfloor\sqrt{N-a^2}\rfloor$. The inner loop over $b=1,\dots,$bmax skips the pairs with std::gcd(a, b) != 1, so that only primitive pairs remain, and adds 2 * a times sigma_prefix(limit / norm, …) with norm $=a^2+b^2$; integer division of positive numbers rounds down, so limit / norm is $\lfloor N/m\rfloor$.

sigma_prefix(n, memo) returns a stored value if memo already holds $T(n)$. Otherwise it runs the blocks of Section 6: with q = n / l and r = n / q it adds q times the interval sum (l + r) * (r - l + 1) / 2, which is exact because $(l+r)+(r-l+1)=2r+1$ is odd, so one of the two factors is even, and it continues with l = r + 1; at the end it stores $T(n)$ in memo. The product of the two factors is at most $2M^2$, and for $M=10^8$ it stays below $10^{16}$; by Section 8 the total fits in 64 bits as well.

The Python program has the same two functions, with a dictionary as the table, math.isqrt for the exact integer square root and math.gcd; its integers have unlimited precision, and it runs the fixed limit $10^8$ without options or checks. The Java program also fixes the limit and has no checks; it uses long sums, its own Euclidean gcd loop and Math.sqrt, and its sigmaPrefix has no table, so it evaluates $T(\lfloor N/m\rfloor)$ afresh for each of the $47\,746\,495$ primitive pairs.

Complexity and Verification

The double loop visits every pair $a,b\ge1$ with $a^2+b^2\le N$. Each such pair owns the unit square with upper right corner $(a,b)$, and these squares lie in the quarter disk of radius $\sqrt N$, so there are at most $\pi N/4$ pairs; for $N=10^8$ there are $78\,529\,763$, against $\pi N/4\approx78\,539\,816$, close to the area as in Gauss's circle problem. ^[8] Each pair costs one gcd. Euclid's algorithm, which the Java gcd spells out, needs at most $20$ division steps here by Lamé's theorem, because the smaller of $a$ and $b$ is at most $7071$ and therefore has at most four decimal digits. ^[9] About $6/\pi^2$ of the pairs are coprime, as the primitive circle problem predicts: $47\,746\,495$ here, against $3N/(2\pi)\approx47\,746\,483$. Counting arithmetic on these 64-bit-sized numbers and each table operation as constant time, the running time is therefore $O(N\log N)$ in all three languages.

With the table, each different argument is evaluated once. An argument $\lfloor N/m\rfloor$ with $m\le\sqrt N$, including $N$ itself as the case $m=1$, needs at most $2\sqrt{N/m}$ blocks, which sum to at most $4N^{3/4}$ over all such $m$, and the fewer than $\sqrt N$ arguments below $\sqrt N$ need fewer than $2N^{1/4}$ blocks each, so all evaluations together cost $O(N^{3/4})$. For $N=10^8$ they take $1\,330\,082$ block steps, negligible beside the gcd work, and the table holds $7\,438$ entries. By Section 7 it never holds more than $2\sqrt N$, so C++ and Python need $O(\sqrt N)$ memory, and Java, which keeps no table, $O(1)$. Without the table, as in Java, every primitive pair costs up to $2\sqrt{N/m}$ blocks. On the unit square of $(a,b)$ the distance $r$ from $0$ is at most $\sqrt m$, so $1/\sqrt m$ is at most the integral of $1/r$ over that square, and the integral of $1/r$ over the quarter disk is $(\pi/2)\sqrt N$. Java therefore makes at most $\pi N+2\sqrt N$ block steps, $T(N)$ included. The measured $126\,286\,416$ are of the same order as the loop itself, so the table saves a constant factor, not an order of growth. On a laptop with an Apple M4 Max, the C++ program took $1.3$ s including its checks, the Java program $2.2$ s and the Python program $7.9$ s.

The C++ checks test $S(5)=35$ and $S(10^5)=17924657155$. For this article, a direct search over the definition, which tests every $x+yi$ with $x\gt0$ and $x^2+y^2\le n^2$ by the criterion of Section 1, reproduced the whole table of the statement, and the formula of Section 4 agreed with it for every $N\le400$. The claim of Section 3 was tested for every $n\le400$ and every primitive pair and multiplier with $d^2(a^2+b^2)\le n^2$. The values $s(1),\dots,s(30)$ agree with OEIS A078930, whose comments describe it as the sum of the divisors with positive real part. ^[10] The C++ program prints $1$, $6$, $10$, $23$, $35$ and $161$ for the limits $1$ to $5$ and $10$, as the direct search does, and the three programs print the same number for $N=10^8$.

Footnotes and References

The references below supply the background facts used above. The description of the divisors by primitive pairs, the reduction to the sums $T(\lfloor N/m\rfloor)$ and the operation counts are derived explicitly in this article.

Project Euler 153 — Investigating Gaussian Integers. The official statement defines the Gaussian integers and the Gaussian divisors of a rational integer n, lists the divisors with positive real part and their sums s(n) for n = 1 to 5, gives the sums 35 for n ≤ 5 and 17924657155 for n ≤ 10⁵, and asks for the sum of s(n) over n ≤ 10⁸.
Wikipedia — Gaussian integer. The Gaussian integers are the complex numbers a + bi with integers a and b; the conjugate of a + bi is a − bi, and the norm N(a + bi) = a² + b², the product of a number and its conjugate, is multiplicative. The units are 1, −1, i and −i.
Wikipedia — Bézout's identity. For integers a and b with greatest common divisor d there are integers x and y with ax + by = d; for coprime a and b this is ax + by = 1, the form used in Section 3.
OEIS A024916. The sequence a(n) = Σ k⌊n/k⌋ over k = 1, …, n, which is also the sum of σ(k) over k ≤ n. Its terms begin 1, 4, 8, 15, 21, 33, 41, 56, 69, 87, and the entry gives the asymptotic formula a(n) = π²n²/12 + O(n log n).
Wikipedia — Floor and ceiling functions. The section on nested divisions states ⌊⌊x/m⌋/n⌋ = ⌊x/(mn)⌋ for a positive integer n; Section 4 proves the case needed here by counting multiples.
Wikipedia — Divisor summatory function. D(x) = Σ⌊x/k⌋ over k ≤ x counts the lattice points under the hyperbola jk = x, and Dirichlet's formula D(x) = 2Σ⌊x/k⌋ − u², summed over k ≤ u with u = ⌊√x⌋, computes it in O(√x) time.
C++ — Euler153.cpp. The linked immutable C++ revision contains sigma_prefix, solve and the two checks in run_checkpoints. It is the implementation record for this article.
Wikipedia — Gauss circle problem. The number of integer points (m, n) with m² + n² ≤ r² is πr² + E(r), and Gauss showed that the error E(r) is at most a constant times r. The primitive circle problem counts the coprime solutions, V(r) = (6/π)r² + O(r^(1+ε)), in line with the probability 6/π² that two integers are coprime.
Wikipedia — Euclidean algorithm. Lamé's theorem of 1844: the algorithm never needs more steps than five times the number of decimal digits of the smaller number.
OEIS A078930. The sequence is defined through the sum r + si of the distinct first-quadrant Gaussian integers dividing n, as the values r + s; a comment states that it also equals the sum of all divisors of n in ℤ[i] with positive real part, and the entry links Project Euler 153. Its terms begin 1, 5, 4, 13, 12, 20, 8, 29, 13, 56.

Problem 153 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <unordered_map>

namespace {

using u64 = std::uint64_t;

struct Options {
    int limit = 100000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 1;
}

u64 sigma_prefix(const int n, std::unordered_map<int, u64>& memo) {
    const auto it = memo.find(n);
    if (it != memo.end()) {
        return it->second;
    }

    u64 sum = 0;
    int l = 1;
    while (l <= n) {
        const int q = n / l;
        const int r = n / q;
        const u64 interval_sum = (static_cast<u64>(l) + static_cast<u64>(r)) *
                                 static_cast<u64>(r - l + 1) / 2ULL;
        sum += static_cast<u64>(q) * interval_sum;
        l = r + 1;
    }
    memo.emplace(n, sum);
    return sum;
}

u64 solve(const int limit) {
    std::unordered_map<int, u64> sigma_cache;
    sigma_cache.reserve(32768);

    u64 answer = sigma_prefix(limit, sigma_cache);
    for (int a = 1; static_cast<u64>(a) * static_cast<u64>(a) <= static_cast<u64>(limit); ++a) {
        const int aa = a * a;
        const int bmax = static_cast<int>(std::sqrt(static_cast<double>(limit - aa)));
        for (int b = 1; b <= bmax; ++b) {
            if (std::gcd(a, b) != 1) {
                continue;
            }
            const int norm = aa + b * b;
            answer += static_cast<u64>(2 * a) * sigma_prefix(limit / norm, sigma_cache);
        }
    }
    return answer;
}

bool run_checkpoints() {
    if (solve(5) != 35ULL) {
        std::cerr << "Checkpoint failed for limit 5" << '\n';
        return false;
    }
    if (solve(100000) != 17924657155ULL) {
        std::cerr << "Checkpoint failed for limit 100000" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

import math

def sigma_prefix(n, memo):
    if n in memo:
        return memo[n]
    s = 0
    l = 1
    while l <= n:
        q = n // l
        r = n // q
        interval_sum = (l + r) * (r - l + 1) // 2
        s += q * interval_sum
        l = r + 1
    memo[n] = s
    return s

def solve():
    limit = 100000000
    memo = {}
    answer = sigma_prefix(limit, memo)
    a = 1
    while a * a <= limit:
        aa = a * a
        bmax = math.isqrt(limit - aa)
        for b in range(1, bmax + 1):
            if math.gcd(a, b) != 1:
                continue
            norm = aa + b * b
            answer += 2 * a * sigma_prefix(limit // norm, memo)
        a += 1
    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

public class Euler153 {
    static long sigmaPrefix(int n) {
        long s = 0;
        int l = 1;
        while (l <= n) {
            int q = n / l, r = n / q;
            long intSum = ((long) l + r) * (r - l + 1) / 2;
            s += (long) q * intSum;
            l = r + 1;
        }
        return s;
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    public static void main(String[] args) {
        int limit = 100000000;
        long answer = sigmaPrefix(limit);
        for (int a = 1; (long) a * a <= limit; a++) {
            int bmax = (int) Math.sqrt(limit - (long) a * a);
            for (int b = 1; b <= bmax; b++) {
                if (gcd(a, b) != 1)
                    continue;
                int norm = a * a + b * b;
                answer += 2L * a * sigmaPrefix(limit / norm);
            }
        }
        System.out.println(answer);
    }
}