Problem 401: Sum of Squares of Divisors

Project Euler Problem 401 Solution

EulerSolve provides an optimized solution for Project Euler Problem 401, Sum of Squares of Divisors, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We need the summatory function of squared divisor sums, $$\Sigma_2(N)=\sum_{n=1}^{N}\sigma_2(n),\qquad \sigma_2(n)=\sum_{d\mid n} d^2,$$ for the very large input $N=10^{15}$, and we only need the result modulo $10^9$. A direct evaluation of $\sigma_2(n)$ for every $n\le N$ is far too slow, so the implementation reorganizes the summation and then skips large ranges at once. Mathematical Approach Step 1: Exchange the Order of Summation Start from the definition of $\sigma_2(n)$. If we sum over all $n\le N$, every divisor $d$ contributes its square $d^2$ once for each multiple of $d$ up to $N$. Therefore $$\Sigma_2(N)=\sum_{n=1}^{N}\sum_{d\mid n} d^2=\sum_{d=1}^{N} d^2\left\lfloor\frac{N}{d}\right\rfloor.$$ This identity is the key simplification. Instead of asking for the divisors of each $n$, we ask how many integers up to $N$ are divisible by a fixed $d$. That count is exactly $\left\lfloor N/d\right\rfloor$. Step 2: Group Equal Quotients into Blocks The transformed sum still runs over $d=1,2,\dots,N$, but the value $$q(d)=\left\lfloor\frac{N}{d}\right\rfloor$$ does not change at every step....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, let $\sigma_2(n)$ be the sum of the squares of the divisors of $n$ and let $\Sigma_2$ be its summatory function, written $\operatorname{SIGMA}_2$ in the statement:

$$\begin{aligned}\sigma_2(n)&=\sum_{d\mid n}d^2,\\ \Sigma_2(N)&=\sum_{n=1}^{N}\sigma_2(n).\end{aligned}$$

The divisors of $6$ are $1,2,3$ and $6$, so $\sigma_2(6)=1+4+9+36=50$. The first six values of $\Sigma_2$ are $1,6,16,37,63$ and $113$. We must find $\Sigma_2(10^{15})$ modulo $10^9$. ^[1] ^[2]

The whole solution in one sentence: count how often each divisor square is used, gather the divisors that share a count into about $2\sqrt N$ blocks, and add every block in one step with the closed formula for a sum of consecutive squares.

Mathematical Approach

1. The two functions on small numbers

The divisors, then $\sigma_2$, then the running total:

$n=1$: $1$; $\sigma_2=1$; $\Sigma_2=1$.
$n=2$: $1,2$; $\sigma_2=5$; $\Sigma_2=6$.
$n=3$: $1,3$; $\sigma_2=10$; $\Sigma_2=16$.
$n=4$: $1,2,4$; $\sigma_2=21$; $\Sigma_2=37$.
$n=5$: $1,5$; $\sigma_2=26$; $\Sigma_2=63$.
$n=6$: $1,2,3,6$; $\sigma_2=50$; $\Sigma_2=113$.
$n=7$: $1,7$; $\sigma_2=50$; $\Sigma_2=163$.
$n=8$: $1,2,4,8$; $\sigma_2=85$; $\Sigma_2=248$.
$n=9$: $1,3,9$; $\sigma_2=91$; $\Sigma_2=339$.
$n=10$: $1,2,5,10$; $\sigma_2=130$; $\Sigma_2=469$.

The first six totals are the statement's, confirming that the sum starts at $n=1$ and includes $n=N$. The function $\sigma_2$ is irregular — $\sigma_2(6)=\sigma_2(7)=50$, although $6$ has four divisors and $7$ only two — so the sum must be attacked as a whole.

2. Counting each divisor by its multiples

Write $\Sigma_2(N)$ as a sum over pairs $(d,n)$, one for each $n\le N$ and each divisor $d$ of $n$, with weight $d^2$. Summing over $n$ first gives the definition; summing over $d$ first asks how many $n\le N$ are multiples of a fixed $d$, and of $d,2d,3d,\dots$ exactly $\lfloor N/d\rfloor$ do not exceed $N$. Counting both ways proves

$$\Sigma_2(N)=\sum_{d=1}^{N}d^2\left\lfloor\frac{N}{d}\right\rfloor.$$

The Dirichlet convolution $f*g$ sums $f(d)g(e)$ over the pairs with $de=n$, so with $\mathrm{Id}_2(d)=d^2$ and the constant $1$ this says $\sigma_2=\mathrm{Id}_2*1$; over $n\le N$ it collects the pairs with $de\le N$. ^[3] For $N=6$ the new form reads

$$1\cdot6+4\cdot3+9\cdot2+16\cdot1+25\cdot1+36\cdot1,$$

that is $6+12+18+16+25+36=113$, the value from the table. No factorisation is needed any more, only one division per term — but there are still $N$ terms, and $10^{15}$ of them is out of reach.

3. Sums of consecutive squares in closed form

The next step needs partial sums of squares:

$$S_2(x)=\sum_{k=1}^{x}k^2=\frac{x(x+1)(2x+1)}{6},$$

with $S_2(0)=0$. ^[4] A short proof: $(k+1)^3-k^3=3k^2+3k+1$ holds for every $k$, and summing it for $k=1,\dots,x$ telescopes the left side to $(x+1)^3-1$, so

$$(x+1)^3-1=3S_2(x)+3\cdot\frac{x(x+1)}{2}+x,$$

and solving for $S_2(x)$ gives the formula; for instance $1+4+9+16+25+36=91$ and $6\cdot7\cdot13/6=91$. Any run of consecutive squares is a difference of two such values:

$$\sum_{d=l}^{r}d^2=S_2(r)-S_2(l-1).$$

4. Blocks on which the quotient is constant

For $N=10$ the quotients $\lfloor10/d\rfloor$, $d=1,\dots,10$, are $10,5,3,2,2,1,1,1,1,1$: all different for small $d$, then repeating in long blocks. If the last index of a block can be found directly, the whole block is added in one step.

Claim. Let $1\le l\le N$, put $q=\lfloor N/l\rfloor$ and $r=\lfloor N/q\rfloor$. Then $l\le r\le N$, $\lfloor N/d\rfloor=q$ for every $d$ with $l\le d\le r$, and $\lfloor N/(r+1)\rfloor\lt q$.

Proof. Since $q\le N/l$ we get $l\le N/q$, and $l$ is an integer, so $l\le r$; also $q\ge1$ gives $r\le N$. The map $d\mapsto\lfloor N/d\rfloor$ never increases, so for $l\le d\le r$ it lies between $\lfloor N/r\rfloor$ and $\lfloor N/l\rfloor=q$. From $r\le N/q$ we get $N/r\ge q$, hence $\lfloor N/r\rfloor\ge q$, so both bounds equal $q$. Finally $r=\lfloor N/q\rfloor$ means $r+1\gt N/q$, hence $\lfloor N/(r+1)\rfloor\lt q$. ∎

For an arbitrary $l$, $[l,r]$ is only the tail of its block: at $N=10$, $l=7$ gives $r=10$ while the block of quotient $1$ starts at $6$. In the walk this cannot happen: $l=1$ opens the first block, and every later $l$ is the previous $r+1$, where the quotient has just dropped. By induction every visited $[l,r]$ is a complete block, the blocks tile $\{1,\dots,N\}$, and there is one per distinct value of $\lfloor N/d\rfloor$. Each is added at once:

$$\sum_{d=l}^{r}d^2\lfloor N/d\rfloor=q\bigl(S_2(r)-S_2(l-1)\bigr).$$

If the blocks are $[l_1,r_1],\dots,[l_B,r_B]$ with quotients $q_j=\lfloor N/l_j\rfloor$, then

$$\Sigma_2(N)=\sum_{j=1}^{B}q_j\bigl(S_2(r_j)-S_2(l_j-1)\bigr).$$

5. Why there are only about $2\sqrt N$ blocks

Put $s=\lfloor\sqrt N\rfloor$. At most $s$ blocks start at a $d\le s$, there being only $s$ such starting points. For $d\ge s+1$,

$$\left\lfloor\frac{N}{d}\right\rfloor\le\frac{N}{s+1}\lt s+1,$$

because $(s+1)^2\gt N$. So there the quotient is one of $1,\dots,s$: at most $s$ further blocks. Altogether $B\le2s$, and the exact count follows too. For $d\lt s$ the gap $N/d-N/(d+1)=N/(d(d+1))$ exceeds $1$, as $d(d+1)\le s^2-s\lt N$, so the first $s$ quotients strictly decrease and are distinct. Put $t=\lfloor N/(s+1)\rfloor$. For $d\ge s+1$ the same gap is below $1$, since $(s+1)(s+2)\gt N$, so the floors descend from $t$ to $1$ without skipping a value: $t$ further blocks (for $N=1$ that range is empty and $t=0$). If $s^2+s\le N$ then $t=s$ and $\lfloor N/s\rfloor\ge s+1\gt t$; otherwise $t=s-1$ and $\lfloor N/s\rfloor\ge s\gt t$. Either way the groups are disjoint, so $B=s+t$: $2s$ in the first case, $2s-1$ in the second, as a check over every $N\le30000$ confirms. At $N=10^{15}$, $s=31622776$ and $s^2+s=999999993568952\le N$, so the walk takes exactly $63245552$ blocks instead of $10^{15}$ terms. Smaller cases: $62$ blocks for $10^3$, $199$ for $10^4$, $1999$ for $10^6$. The same grouping is the standard way to evaluate $D(x)=\sum_{k\le x}\lfloor x/k\rfloor$. ^[5]

6. A worked example: $N=10$

From $l=1$ the claim of Section 4 produces five blocks; for each, $l$, $q$, $r$, the sum of squares and the product:

$l=1$, $q=10$, $r=\lfloor10/10\rfloor=1$; squares $1$; product $10$.
$l=2$, $q=5$, $r=\lfloor10/5\rfloor=2$; squares $4$; product $20$.
$l=3$, $q=3$, $r=\lfloor10/3\rfloor=3$; squares $9$; product $27$.
$l=4$, $q=2$, $r=\lfloor10/2\rfloor=5$; squares $S_2(5)-S_2(3)=55-14=41$; product $82$.
$l=6$, $q=1$, $r=\lfloor10/1\rfloor=10$; squares $S_2(10)-S_2(5)=385-55=330$; product $330$.

The total is $10+20+27+82+330=469$, the value $\Sigma_2(10)$ listed in Section 1. Here $s=3$ and $s^2+s=12\gt10$, so indeed $B=2s-1=5$.

7. The modulus $10^9$ is not prime, and $6$ is not invertible

Only the result modulo $10^9$ is wanted, and $10^9=2^9\cdot5^9$ is not prime. That matters, because the closed form of $S_2$ divides by $6$, and dividing modulo $m$ needs an inverse of the divisor, which exists exactly when it is coprime to $m$. ^[6] Here $\gcd(6,10^9)=2$, so there is none; directly, $6k\equiv1$ would mean $6k-1=10^9j$, odd on the left and even on the right.

Nor may one reduce first and divide afterwards. At $x=794$, $794\cdot795\cdot1589=1003024470$ reduces to $3024470$, which is not even a multiple of $6$: it is $6\cdot504078+2$. Yet $S_2(794)=167170745$ is a perfectly good integer; reduction destroys the divisibility.

Evaluating $x(x+1)(2x+1)$ exactly and dividing before reducing is correct but too large for fixed-width arithmetic: at $x=10^{15}$ the product has $46$ decimal digits, while $2^{64}$ has $20$ and $2^{128}$ only $39$. Arbitrary-precision integers can take that route; the cancellation below avoids it entirely.

8. Cancelling the $6$ before reducing

The cure is to remove the factors $2$ and $3$ from $x$, $x+1$, $2x+1$ while they are small, and reduce only then. Two facts make this possible.

Exactly one of $x$ and $x+1$ is even, since consecutive integers have different parities, so dividing that one by $2$ is exact.
Exactly one of $x$, $x+1$, $2x+1$ is a multiple of $3$, by the residue of $x$: if $x\equiv0$ then $3\mid x$; if $x\equiv1$ then $2x+1\equiv3\equiv0$; if $x\equiv2$ then $x+1\equiv0$. The other two are then not multiples of $3$.

After these two exact divisions the factors multiply to $S_2(x)$ exactly, each is at most $2x+1$, and each may be reduced modulo $10^9$ on its own. For $x=9$ the triple $(9,10,19)$ becomes $(9,5,19)$ and then $(3,5,19)$, with $3\cdot5\cdot19=285=S_2(9)$; for $x=10$ it becomes $(5,11,21)$ and then $(5,11,7)$, with $5\cdot11\cdot7=385=S_2(10)$.

One detail deserves attention: the implementations halve first and test divisibility by $3$ on the halved value. That is harmless, since $3$ and $2$ are coprime, so $3\mid y/2$ exactly when $3\mid y$; dividing by $3$ likewise preserves parity.

As a side remark, $3$ alone is invertible modulo $10^9$, with $3^{-1}=666666667$, so cancelling only the $2$ and multiplying by that inverse works too, as was checked. The cancellation uses no inverse and is valid for every modulus $m\ge2$; the implementations add limits of their own, discussed below.

9. The formula that is actually evaluated

Collecting Sections 4, 5 and 8, the computation is a single walk over the blocks. The quantities $N$, $l$, $q$ and $r$ stay exact integers, because the floor divisions need them; only the square sums, the block contributions and the running total are reduced modulo $m=10^9$, and $q$ is reduced where it is used as a multiplier. Writing $A_j=S_2(r_j)-S_2(l_j-1)$,

$$\Sigma_2(N)\equiv\sum_{j=1}^{B}q_jA_j\pmod m.$$

Addition and multiplication of residues respect the modulus, so reductions may be applied at every step. Only the subtraction needs care: $S_2(r)-S_2(l-1)$ can come out negative once both terms are reduced, so $m$ is added before the final reduction.

10. Two independent rearrangements, used as checks

The pairs $(d,e)$ with $de\le N$ can be grouped by $e$ instead of by $d$. For fixed $e$ the admissible $d$ are $1,\dots,\lfloor N/e\rfloor$, so

$$\Sigma_2(N)=\sum_{e=1}^{N}S_2\bigl(\lfloor N/e\rfloor\bigr).$$

Section 4 applies to $e$ too: on a block $[l,r]$ with $q=\lfloor N/l\rfloor$ every $e$ contributes the same $S_2(q)$, so the block contributes $(r-l+1)S_2(q)$ — the multiplicity must not be dropped. A third arrangement splits the pairs at $K=\lfloor\sqrt N\rfloor$, with $L=\lfloor N/K\rfloor$. Every pair with $de\le N$ has $d\le K$ or $e\le L$, since $d\ge K+1$, $e\ge L+1$ would give $de\gt(K+1)N/K\gt N$; and the whole rectangle $d\le K$, $e\le L$ is admissible, as $KL\le N$. It is counted twice, with weight $L\,S_2(K)$, so

$$\begin{aligned}\Sigma_2(N)=&\sum_{d=1}^{K}d^2\lfloor N/d\rfloor\\ &+\sum_{e=1}^{L}S_2(\lfloor N/e\rfloor)-L\,S_2(K).\end{aligned}$$

This is the hyperbola split that gives $D(x)=2\sum_{k\le u}\lfloor x/k\rfloor-u^2$ for the plain divisor count. ^[5] All three agree; see the verification section.

11. How large is the true value

The second form explains why the problem asks for a residue. For $t\ge1$ the closed form gives $S_2(\lfloor t\rfloor)=t^3/3+O(t^2)$, because $\lfloor t\rfloor=t+O(1)$. Substituting $t=N/e$ and summing over $e\le N$, the errors total $O(N^2)$ since $\sum_{e\ge1}e^{-2}$ converges, and completing the main sum to all $e$ costs only $O(N)$ since $\sum_{e\gt N}e^{-3}\le1/(2N^2)$. Hence

$$\Sigma_2(N)=\frac{\zeta(3)}{3}N^3+O(N^2),$$

with $\zeta(3)=1.2020569\ldots$, Apéry's constant. ^[7] The leading term is already close at modest sizes: $\Sigma_2(10^6)=400686363385965077$ against $\zeta(3)\cdot10^{18}/3=4.006856\ldots\cdot10^{17}$, a relative error of about $1.8\cdot10^{-6}$. At $N=10^{15}$ the value is about $4.0\cdot10^{44}$, a number of $45$ decimal digits — hence the residue.

How the Code Works

The C++ program ^[10] keeps its parameters in an Options record: $n=10^{15}$, $\mathrm{mod}=10^9$ and a flag for the self-checks. parse_u64_after_prefix and parse_arguments read the overrides --n=, --mod= and --skip-checkpoints and reject anything else, $n\lt1$ or $\mathrm{mod}\lt2$. The parser does not detect overflow, and $2n+1$, right + mod and ans + term stay inside 64 bits only while $n\lt2^{63}$ and $2\,\mathrm{mod}\le2^{64}$; with --mod=18446744073709551615 and $n=1$ the program returns $0$ instead of $1$. The defaults are far inside those bounds.

sum_squares_mod(n, mod) is Section 8. It returns $0$ for $n=0$, which makes $S_2(l-1)$ legitimate in the first block, and otherwise sets $a=n$, $b=n+1$, $c=2n+1$, halves whichever of $a$, $b$ is even, then divides by $3$ whichever of $a$, $b$, $c$ is a multiple of it, testing $a$, then $b$, falling back to $c$. The factors are reduced modulo mod and multiplied in two steps through __uint128_t, a GCC and Clang extension for 128-bit integers. ^[8] sum_squares_range_mod(l, r, mod) is Section 3 with the guarded subtraction of Section 9, returning (right + mod - left) % mod.

solve(n, mod) is the block walk of Sections 4 and 9: from l = 1 it computes q = n / l and r = n / q, takes the block's sum of squares, multiplies by q % mod, adds it to ans and continues with l = r + 1. It uses no array and enumerates no divisors; only the checkpoint routine below does.

brute_sigma2_sum(n) is the slow reference: for each $x\le n$ it tries every $d$ with $d^2\le x$ and on a hit adds $d^2$ and $(x/d)^2$, skipping the second when $d=x/d$, in an exact sum without any modulus. run_checkpoints compares solve with the six statement values and with brute_sigma2_sum(1000) reduced modulo $10^9$; main runs it unless switched off, exits with status $2$ on a failure, and otherwise prints solve.

The Python version packs the same pieces into nested functions sum_sq_mod and sum_sq_range inside solve, with $N$ and MOD fixed in the body. Its integers have unlimited precision, ^[9] so it could have formed $x(x+1)(2x+1)$ and divided by $6$; it mirrors the cancellation instead. With a positive modulus its remainder lies between $0$ and $m-1$, so the subtraction needs no repair.

The Java version uses plain long throughout: sumSquaresMod, sumSquaresRangeMod with an explicit if (diff < 0) diff += mod;, and the same loop in solve. No wide type is needed: the largest intermediate is a product of two reduced residues, at most $(10^9-1)^2=999999998000000001$, well below $2^{63}-1=9223372036854775807$. ^[9] Only the C++ file carries the checkpoints.

Complexity and Verification

The walk visits one block per distinct value of $\lfloor N/d\rfloor$, so by Section 5 it runs in $\Theta(\sqrt N)$ time: exactly $63245552$ iterations at $N=10^{15}$, each a constant number of integer divisions and modular multiplications. Memory is $O(1)$: the loop variables and a few residues, no sieve, no table. The slow reference is $O(n^{3/2})$ and runs only at $n=1000$. Measured here, C++ needs about $1.8$ seconds, Java about $0.6$ and Python about $37$, the last because $63$ million interpreted iterations dominate everything else. All three print the same value.

The packaged checks are the six statement values $\Sigma_2(1)=1$ through $\Sigma_2(6)=113$ and the comparison with the brute force at $n=1000$. That comparison is sharper than it looks: $\Sigma_2(1000)=401382971$ is below $10^9$, so the modulus does not act and the two sides agree as integers. The accumulator of brute_sigma2_sum is an unsigned 64-bit integer, and by Section 11 the exact sum passes $2^{64}$ near $n\approx3.58\cdot10^6$, so the reference cannot be moved to much larger $n$ unchanged.

Independently of the published files, the following was checked for this article. The block walk matches the plain sum $\sum_{d\le N}d^2\lfloor N/d\rfloor$ at $N=10^5$, $10^6$ and $10^7$. The block walk, the grouping by $e$, the hyperbola split and a direct divisor-by-divisor brute force give the same residue for every $N\le3000$, and at $N=10^{15}$ the three fast arrangements agree, taking about $0.47$, $0.39$ and $0.20$ seconds in a separate compiled check. The cancellation of Section 8 matched the exact $x(x+1)(2x+1)/6$ modulo $10^9$ for every $x\lt30000$ and for $20004$ random and boundary values up to $10^{15}$, as did the variant using $3^{-1}=666666667$. The block claim was verified for all $N\le2000$ and the exact block count for all $N\le30000$.

Footnotes and References

The references supply the background facts used above; the claims of Sections 4, 5, 8 and 10 are derived in this article.

Project Euler 401 — Sum of Squares of Divisors. The official statement defines sigma_2(n) as the sum of the squares of the divisors of n, notes that the divisors of 6 give 50, introduces SIGMA_2 as its summatory function, lists its first six values as 1, 6, 16, 37, 63 and 113, and asks for SIGMA_2(10¹⁵) modulo 10⁹.
OEIS A064602 — partial sums of sigma_2. The entry is the sum of sigma_2(j) for j = 1 to n, with sigma_2 = A001157, and its first terms are 1, 6, 16, 37, 63, 113, 163, 248, 339, 469. These are the values tabulated in Section 1.
Wikipedia — Dirichlet convolution. The convolution of f and g is the sum of f(a)g(b) over all pairs of positive integers with ab = n, and σ_k = Id_k ∗ 1 is among the article's examples, as used in Section 2.
Wikipedia — Square pyramidal number. The article gives the sum of the first n squares as n(n + 1)(2n + 1)/6, a special case of Faulhaber's formula provable by induction, and records that Archimedes used the sum as a lemma. Section 3 proves it by telescoping.
Wikipedia — Divisor summatory function. The article writes D(x) as the number of pairs (j, k) with jk ≤ x, gives the grouped form D(x) = 2∑_{k≤u}⌊x/k⌋ − u², u = ⌊√x⌋, as a way to compute it in O(√x) time, and attributes its O(√x) error term to the Dirichlet hyperbola method. Sections 5 and 10 are the σ₂ analogues.
Wikipedia — Modular multiplicative inverse. The article states that an inverse of a modulo m exists if and only if gcd(a, m) = 1, that is, exactly when a and m are relatively prime. Section 7 applies it to a = 6 and m = 10⁹, whose gcd is 2.
Wikipedia — Apéry's constant. It is ζ(3), the value of the Riemann zeta function at 3 and the sum of the reciprocals of the cubes, approximately 1.2020569031595942854; Apéry proved in 1978 that it is irrational. Section 11 uses the value.
GCC manual — 128-bit integers. The page documents __int128 and unsigned __int128 as an extension for targets with an integer mode wide enough to hold 128 bits; the C++ version uses the unsigned form for its modular products.
Python — numeric types states that integers have unlimited precision, so the Python version could divide the full product by 6 directly. Java — Long documents MAX_VALUE as 2⁶³ − 1 for the 64-bit signed long, the bound used in the code section.
C++ — Euler401.cpp. The linked immutable revision contains sum_squares_mod, sum_squares_range_mod, solve, brute_sigma2_sum and run_checkpoints, and is the implementation record for this article; the extra comparisons of the verification section are not part of it.

Problem 401 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;

struct Options {
    u64 n = 1'000'000'000'000'000ULL;
    u64 mod = 1'000'000'000ULL;
    bool run_checkpoints = true;
};

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

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

u64 sum_squares_mod(u64 n, const u64 mod) {
    if (n == 0ULL) {
        return 0ULL;
    }
    u64 a = n;
    u64 b = n + 1ULL;
    u64 c = 2ULL * n + 1ULL;

    if ((a & 1ULL) == 0ULL) {
        a /= 2ULL;
    } else {
        b /= 2ULL;
    }
    if (a % 3ULL == 0ULL) {
        a /= 3ULL;
    } else if (b % 3ULL == 0ULL) {
        b /= 3ULL;
    } else {
        c /= 3ULL;
    }

    const __uint128_t x = static_cast<__uint128_t>(a % mod) * static_cast<__uint128_t>(b % mod);
    const __uint128_t y = (x % mod) * static_cast<__uint128_t>(c % mod);
    return static_cast<u64>(y % mod);
}

u64 sum_squares_range_mod(const u64 l, const u64 r, const u64 mod) {
    const u64 right = sum_squares_mod(r, mod);
    const u64 left = sum_squares_mod(l - 1ULL, mod);
    return (right + mod - left) % mod;
}

u64 solve(const u64 n, const u64 mod) {
    u64 ans = 0ULL;
    for (u64 l = 1ULL; l <= n;) {
        const u64 q = n / l;
        const u64 r = n / q;
        const u64 sum_d2 = sum_squares_range_mod(l, r, mod);
        const __uint128_t term = static_cast<__uint128_t>(q % mod) * static_cast<__uint128_t>(sum_d2);
        ans = (ans + static_cast<u64>(term % mod)) % mod;
        l = r + 1ULL;
    }
    return ans;
}

u64 brute_sigma2_sum(const u64 n) {
    u64 s = 0ULL;
    for (u64 x = 1ULL; x <= n; ++x) {
        u64 sigma2 = 0ULL;
        for (u64 d = 1ULL; d * d <= x; ++d) {
            if (x % d != 0ULL) {
                continue;
            }
            sigma2 += d * d;
            const u64 e = x / d;
            if (e != d) {
                sigma2 += e * e;
            }
        }
        s += sigma2;
    }
    return s;
}

bool run_checkpoints() {
    const u64 mod = 1'000'000'000ULL;
    if (solve(1ULL, mod) != 1ULL ||
        solve(2ULL, mod) != 6ULL ||
        solve(3ULL, mod) != 16ULL ||
        solve(4ULL, mod) != 37ULL ||
        solve(5ULL, mod) != 63ULL ||
        solve(6ULL, mod) != 113ULL) {
        std::cerr << "Checkpoint failed: first SIGMA2 values" << '\n';
        return false;
    }
    if (solve(1000ULL, mod) != brute_sigma2_sum(1000ULL) % mod) {
        std::cerr << "Checkpoint failed: brute-force cross-check for n=1000" << '\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.n, options.mod) << '\n';
    return 0;
}

Python

def solve():
    N = 10**15
    MOD = 10**9

    def sum_sq_mod(n, mod):
        if n == 0: return 0
        a, b, c = n, n+1, 2*n+1
        if a % 2 == 0: a //= 2
        else: b //= 2
        if a % 3 == 0: a //= 3
        elif b % 3 == 0: b //= 3
        else: c //= 3
        return (a % mod) * (b % mod) % mod * (c % mod) % mod

    def sum_sq_range(l, r, mod):
        return (sum_sq_mod(r, mod) - sum_sq_mod(l-1, mod)) % mod

    ans = 0
    l = 1
    while l <= N:
        q = N // l
        r = N // q
        term = (q % MOD) * sum_sq_range(l, r, MOD) % MOD
        ans = (ans + term) % MOD
        l = r + 1

    return str(ans)

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

Java

public class Euler401 {
    static long sumSquaresMod(long n, long mod) {
        if (n == 0)
            return 0;
        long a = n;
        long b = n + 1;
        long c = 2 * n + 1;

        if ((a & 1) == 0)
            a /= 2;
        else
            b /= 2;

        if (a % 3 == 0)
            a /= 3;
        else if (b % 3 == 0)
            b /= 3;
        else
            c /= 3;

        long x = ((a % mod) * (b % mod)) % mod;
        long y = (x * (c % mod)) % mod;
        return y;
    }

    static long sumSquaresRangeMod(long l, long r, long mod) {
        long right = sumSquaresMod(r, mod);
        long left = sumSquaresMod(l - 1, mod);
        long diff = right - left;
        if (diff < 0)
            diff += mod;
        return diff % mod;
    }

    static String solve() {
        long n = 1000000000000000L;
        long mod = 1000000000L;
        long ans = 0;
        long l = 1;

        while (l <= n) {
            long q = n / l;
            long r = n / q;
            long sumD2 = sumSquaresRangeMod(l, r, mod);
            long term = ((q % mod) * sumD2) % mod;
            ans = (ans + term) % mod;
            l = r + 1;
        }

        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}