Detailed mathematical approach

Problem Summary

The divisor-based quantity in this problem can be reindexed as an ordered factor-pair sum:

$S(N)=\sum_{n=1}^{N}\sum_{d\mid n}\gcd\!\left(d,\frac{n}{d}\right)=\sum_{ab\le N}\gcd(a,b).$

So the task is to add \(\gcd(a,b)\) over all ordered pairs \((a,b)\) whose product is at most \(N\). At \(N=10^{15}\), direct enumeration is impossible, so the solution rewrites the sum using Euler's totient function and the divisor summatory function.

Mathematical Approach

We work with the equivalent form

$S(N)=\sum_{ab\le N}\gcd(a,b).$

The efficient formula comes from separating the common-divisor contribution, then grouping repeated floor values so the large range can be handled in blocks.

Step 1: Reindex the Divisor Pairs

For each \(n\) and each divisor \(d\mid n\), the pair \((d,n/d)\) is an ordered factor pair whose product is exactly \(n\). Summing over all \(n\le N\) and all divisors therefore gives exactly the same set of terms as summing over all ordered pairs \((a,b)\) with \(ab\le N\).

This is the first important simplification: the original divisor language becomes a two-variable summation problem with a clean product bound.

Step 2: Expand \(\gcd\) with Euler's Totient Function

Use the classical identity

$\sum_{t\mid m}\varphi(t)=m.$

Applying it with \(m=\gcd(a,b)\) gives

$\gcd(a,b)=\sum_{t\mid \gcd(a,b)}\varphi(t).$

Substitute this into the sum and swap the order of summation:

$S(N)=\sum_{ab\le N}\sum_{t\mid a,\ t\mid b}\varphi(t).$

Now write \(a=tx\) and \(b=ty\). The condition becomes

$t^2xy\le N.$

Hence

$S(N)=\sum_{t\le \sqrt N}\varphi(t)\,D\!\left(\left\lfloor\frac{N}{t^2}\right\rfloor\right),$

where \(D(m)\) counts ordered pairs \((x,y)\) with \(xy\le m\).

Step 3: Introduce the Divisor Summatory Function

The pair-counting function is

$D(m)=\#\{(x,y)\in\mathbb{Z}_{>0}^2:xy\le m\}.$

For a fixed \(x\), there are exactly \(\left\lfloor m/x\right\rfloor\) valid choices for \(y\). Therefore

$D(m)=\sum_{x=1}^{m}\left\lfloor\frac{m}{x}\right\rfloor.$

Equivalently, every integer \(n\le m\) contributes once for each divisor pair \((x,y)\) with \(xy=n\), so

$D(m)=\sum_{n=1}^{m}\tau(n),$

where \(\tau(n)\) is the divisor-count function. The whole problem is now a single weighted sum of totients times divisor-summatory values.

Step 4: Evaluate \(D(m)\) by Quotient Blocks

A naive loop up to \(m\) is still too slow when \(m\) is large. The key observation is that the quotient

$q=\left\lfloor\frac{m}{\ell}\right\rfloor$

stays constant on an interval \([\ell,r]\), where

$r=\left\lfloor\frac{m}{q}\right\rfloor.$

So one whole block contributes

$q\,(r-\ell+1).$

Advancing from one block to the next yields

$D(m)=\sum_{\text{blocks }[\ell,r]}\left\lfloor\frac{m}{\ell}\right\rfloor(r-\ell+1),$

and the number of blocks is only \(O(\sqrt m)\). This is the hyperbola-style speedup used whenever \(D(m)\) has to be evaluated for a large argument.

Step 5: Split the Outer Sum at \(N^{1/3}\)

Let

$K=\left\lfloor N^{1/3}\right\rfloor,\qquad M=\left\lfloor\frac{N}{K^2}\right\rfloor.$

Then the outer sum naturally splits into two regimes.

For \(t>K\), we have \(\left\lfloor N/t^2\right\rfloor\le M\), so all needed \(D(m)\) values are small. The implementation precomputes \(\tau(1),\dots,\tau(M)\), takes prefix sums, and obtains every small \(D(m)\) by table lookup.

For \(t\le K\), the argument \(\left\lfloor N/t^2\right\rfloor\) is large, but many consecutive \(t\) share the same value. If

$v=\left\lfloor\frac{N}{t^2}\right\rfloor,$

then the largest index with the same quotient is

$r=\left\lfloor\sqrt{\frac{N}{v}}\right\rfloor.$

That whole interval contributes

$D(v)\sum_{u=t}^{r}\varphi(u).$

Prefix sums of \(\varphi\) make the interval sum an \(O(1)\) lookup, so the expensive work is concentrated only in the distinct large \(D(v)\) calls.

Worked Example: \(N=10\)

Here \(\lfloor\sqrt{10}\rfloor=3\), so the transformed formula becomes

$S(10)=\sum_{t=1}^{3}\varphi(t)\,D\!\left(\left\lfloor\frac{10}{t^2}\right\rfloor\right).$

The required totients are

$\varphi(1)=1,\qquad \varphi(2)=1,\qquad \varphi(3)=2.$

The required divisor-summatory values are

$D(10)=\sum_{x=1}^{10}\left\lfloor\frac{10}{x}\right\rfloor=27,\qquad D(2)=3,\qquad D(1)=1.$

Therefore

$S(10)=1\cdot 27+1\cdot 3+2\cdot 1=32,$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same number-theoretic formula. They first compute the integer square-root limit \(\lfloor\sqrt N\rfloor\) and the cube-root split point \(\lfloor N^{1/3}\rfloor\). Next they build a totient sieve up to \(\lfloor\sqrt N\rfloor\) and store prefix sums so that \(\sum_{u=l}^{r}\varphi(u)\) can be recovered instantly.

For the small-argument region, the implementation precomputes divisor counts up to \(M=\lfloor N/K^2\rfloor\) by visiting multiples of each divisor, then turns those counts into prefix sums to obtain \(D(1),D(2),\dots,D(M)\). For the large-argument region, it groups equal values of \(\left\lfloor N/t^2\right\rfloor\), evaluates each large \(D(v)\) with quotient blocks, and multiplies by the corresponding totient-interval sum.

The C++ and Java implementations parallelize those independent large-group contributions before adding the small tail. The Python implementation is a thin wrapper: it builds the compiled solver when needed, runs it, and returns the printed result. The compiled version also checks the two checkpoints \(S(10)=32\) and \(S(1000)=12776\) before evaluating the full target.

Complexity Analysis

Let \(L=\lfloor\sqrt N\rfloor\), \(K=\lfloor N^{1/3}\rfloor\), and \(M=\lfloor N/K^2\rfloor\). The totient sieve costs \(O(L\log\log L)\) time and \(O(L)\) memory. The small-table preparation for divisor counts costs

$\sum_{d=1}^{M}\left\lfloor\frac{M}{d}\right\rfloor=O(M\log M)$

time and \(O(M)\) memory.

In the large-value region there are only \(O(K)\) distinct quotient groups, and one call to the block method for \(D(v)\) costs \(O(\sqrt v)\). Summing those costs over \(t\le K\) gives

$O\!\left(\sum_{t=1}^{K}\sqrt{\frac{N}{t^2}}\right)=O(\sqrt N\log K).$

So the overall running time is \(O(\sqrt N\log N)\) in the worst case, while memory is dominated by the totient arrays and remains \(O(\sqrt N)\). Parallel execution improves wall-clock time for the grouped large-value part but does not change the asymptotic bound.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=530
2. Euler's totient function and the identity \(\sum_{d\mid n}\varphi(d)=n\): Wikipedia — Euler's totient function
3. Divisor summatory function: Wikipedia — Divisor summatory function
4. Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method
5. Divisor function \(\tau(n)\): Wikipedia — Divisor function