Detailed mathematical approach

Problem Summary

Let \(\sigma(n)\) denote the sum of all positive divisors of \(n\). The problem asks for

$S(N)=\sum_{x=1}^{N}\sum_{y=1}^{N}\sigma(xy)$

with \(N=10^{11}\), reduced modulo \(10^9\). A direct double loop is hopeless, so the solution rewrites \(\sigma(xy)\) into a form that separates the common part of \(x\) and \(y\), then evaluates the resulting sums by quotient blocks.

Mathematical Approach

Step 1: A Möbius Identity for \(\sigma(xy)\)

The key identity used by the implementation is

$\sigma(xy)=\sum_{d\mid \gcd(x,y)} \mu(d)\,d\,\sigma\!\left(\frac{x}{d}\right)\sigma\!\left(\frac{y}{d}\right),$

where \(\mu\) is the Möbius function. Both sides are multiplicative in the pair \((x,y)\), so it is enough to check prime powers. Let \(x=p^a\) and \(y=p^b\).

If \(a=0\) or \(b=0\), then \(\gcd(x,y)=1\) and the identity is immediate. For \(a,b\ge 1\), only \(d=1\) and \(d=p\) contribute, because \(\mu(p^k)=0\) for \(k\ge 2\). Thus

$\begin{aligned} \sum_{d\mid \gcd(p^a,p^b)} \mu(d)\,d\,\sigma\!\left(\frac{p^a}{d}\right)\sigma\!\left(\frac{p^b}{d}\right) &=\sigma(p^a)\sigma(p^b)-p\,\sigma(p^{a-1})\sigma(p^{b-1})\\ &=\frac{(p^{a+1}-1)(p^{b+1}-1)-p(p^a-1)(p^b-1)}{(p-1)^2}\\ &=\frac{p^{a+b+1}-1}{p-1} =\sigma(p^{a+b}). \end{aligned}$

Since \(p^{a+b}=xy\), the identity holds for prime powers and therefore for all \(x,y\).

Step 2: Collapse the Double Sum

Insert the identity into the definition of \(S(N)\):

$S(N)=\sum_{x=1}^{N}\sum_{y=1}^{N}\sum_{d\mid \gcd(x,y)} \mu(d)\,d\,\sigma\!\left(\frac{x}{d}\right)\sigma\!\left(\frac{y}{d}\right).$

Exchange the order of summation and write \(x=da\), \(y=db\). Then \(a,b\le \lfloor N/d\rfloor\), so

$S(N)=\sum_{d=1}^{N}\mu(d)\,d\left(\sum_{a\le N/d}\sigma(a)\right)\left(\sum_{b\le N/d}\sigma(b)\right).$

Define the summatory divisor function

$A(M)=\sum_{n=1}^{M}\sigma(n).$

This gives the single outer sum

$\boxed{S(N)=\sum_{d=1}^{N}\mu(d)\,d\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)^2.}$

Step 3: Fast Evaluation of \(A(M)\)

We start from the divisor definition:

$A(M)=\sum_{n=1}^{M}\sum_{e\mid n} e=\sum_{e=1}^{M} e\left\lfloor\frac{M}{e}\right\rfloor.$

The implementations use an equivalent form that is better suited to quotient blocking. Let

$T(q)=\frac{q(q+1)}{2}.$

Then

$A(M)=\sum_{t=1}^{M} T\!\left(\left\lfloor\frac{M}{t}\right\rfloor\right),$

because for each \(t\), the inner sum \(\sum_{e\le M/t} e\) is exactly \(T(\lfloor M/t\rfloor)\). Now \(\left\lfloor M/t\right\rfloor\) is constant on intervals \(l\le t\le r\) with

$q=\left\lfloor\frac{M}{l}\right\rfloor,\qquad r=\left\lfloor\frac{M}{q}\right\rfloor.$

So one whole block contributes

$\sum_{t=l}^{r} T\!\left(\left\lfloor\frac{M}{t}\right\rfloor\right)=(r-l+1)\,T(q).$

There are only \(O(\sqrt{M})\) such blocks, which is the standard harmonic-sum speedup.

Step 4: Prefix Sums of \(d\,\mu(d)\)

The outer formula needs interval sums of \(\mu(d)\,d\). Define

$W(n)=\sum_{k=1}^{n} k\,\mu(k).$

The implementation computes small values by sieve and large values by memoized recursion. The recurrence comes from the identity

$\sum_{k=1}^{n} k\,W\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right)=1.$

Indeed, after expanding \(W\),

$\sum_{k=1}^{n} k\sum_{m\le n/k} m\,\mu(m)=\sum_{km\le n} km\,\mu(m)=\sum_{t\le n} t\sum_{m\mid t}\mu(m)=1,$

because \(\sum_{m\mid t}\mu(m)=0\) unless \(t=1\). Therefore

$W(n)=1-\sum_{k=2}^{n} k\,W\!\left(\left\lfloor\frac{n}{k}\right\rfloor\right).$

Again, equal quotients are grouped into intervals. For any block \(l\le k\le r\) with \(\lfloor n/k\rfloor=v\), we replace the whole block by

$\left(\sum_{k=l}^{r} k\right)W(v).$

Step 5: Final Harmonic Decomposition

The same block idea is applied to the final sum for \(S(N)\). If \(\lfloor N/d\rfloor=q\) on \(l\le d\le r\), then the entire interval contributes

$A(q)^2\sum_{d=l}^{r} d\,\mu(d)=A(q)^2\bigl(W(r)-W(l-1)\bigr).$

So the computation walks over the \(O(\sqrt N)\) distinct quotient values of \(\lfloor N/d\rfloor\), evaluates \(A(q)\) once per block, and weights it by the corresponding difference of weighted Möbius prefixes.

Worked Example: \(N=3\)

For \(N=3\), the divisor sums are \(\sigma(1)=1\), \(\sigma(2)=3\), \(\sigma(3)=4\), so

$A(3)=1+3+4=8,\qquad A(1)=1.$

The outer formula gives

$S(3)=\mu(1)\cdot 1\cdot A(3)^2+\mu(2)\cdot 2\cdot A(1)^2+\mu(3)\cdot 3\cdot A(1)^2,$

hence

$S(3)=1\cdot 64-2\cdot 1-3\cdot 1=59,$

which matches the small validation target used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they precompute Möbius values up to a practical cutoff and store the prefix sums of \(k\mu(k)\). For larger arguments they reuse the recurrence for \(W(n)\) with memoization, so repeated quotient values are solved only once.

Next they enumerate the quotient blocks of \(\lfloor N/d\rfloor\). For each block they evaluate \(A(q)\) by the triangular-number block formula, square it modulo \(10^9\), multiply by the weighted Möbius difference for that block, and add the result to the running total. The large independent blocks can also be accumulated in parallel.

Complexity Analysis

The direct definition needs \(N^2\) evaluations of \(\sigma(xy)\), which is completely infeasible. After the Möbius decomposition, the outer sum has only \(O(\sqrt N)\) distinct quotient blocks. One evaluation of \(A(M)\) costs \(O(\sqrt M)\), and summing this over the distinct values \(M=\lfloor N/d\rfloor\) gives about \(O(N^{3/4})\) total block work. The weighted Möbius prefix is also handled sublinearly by combining a sieve for small values with memoized quotient recursion for large ones. Memory usage is linear in the chosen sieve cutoff, plus \(O(\sqrt N)\) cached block data.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=439
2. Divisor function \(\sigma(n)\): Wikipedia — Divisor function
3. Möbius function and inversion: Wikipedia — Möbius function
4. Dirichlet convolution: Wikipedia — Dirichlet convolution
5. Harmonic-lemma style floor blocking: cp-algorithms — Number of divisors / sum of divisors