Detailed mathematical approach

Problem Summary

The quantity to evaluate is

$S(N)=\sum_{n \le N} 2^{\Omega(n)},$

where \(\Omega(n)\) is the total number of prime factors of \(n\), counted with multiplicity. For the actual input \(N=10^{14}\), direct enumeration is hopeless, so the solution rewrites the summand as a convolution between the ordinary divisor function and a sparse multiplicative correction supported only on powerful numbers.

Mathematical Approach

The whole method starts from the fact that both \(2^{\Omega(n)}\) and the divisor-counting function \(\tau(n)\) are multiplicative. That makes it natural to search for a second multiplicative function whose convolution with \(\tau\) reproduces the target summand.

Step 1: Introduce a Sparse Multiplicative Function

Define a multiplicative function \(c(n)\) by

$c(1)=1,\qquad c(p)=0,\qquad c(p^e)=2^{e-2}\quad (e\ge 2).$

Now consider the Dirichlet convolution \((\tau * c)(n)\). Since both factors are multiplicative, it is enough to verify the identity on prime powers.

Step 2: Verify the Prime-Power Identity

Let \(n=p^e\). Then

$\begin{aligned} (\tau * c)(p^e) &= \sum_{k=0}^{e} \tau(p^{e-k})\,c(p^k) \\ &= (e+1) + \sum_{k=2}^{e} (e-k+1)\,2^{k-2}. \end{aligned}$

The remaining finite arithmetic-geometric sum simplifies to

$ (e+1) + \sum_{k=2}^{e} (e-k+1)\,2^{k-2} = 2^e. $

Therefore

$ (\tau * c)(p^e)=2^e=2^{\Omega(p^e)}. $

By multiplicativity, this yields the identity

$ 2^{\Omega(n)} = (\tau * c)(n)\quad (n\ge 1). $

Step 3: Turn the Problem into a Summatory Formula

Introduce the divisor summatory function

$D(x)=\sum_{m \le x}\tau(m).$

Using the convolution identity, we get

$\begin{aligned} S(N) &= \sum_{n \le N} (\tau * c)(n) \\ &= \sum_{q \le N} c(q)\sum_{m \le N/q}\tau(m) \\ &= \sum_{q \le N} c(q)\,D\!\left(\left\lfloor \frac{N}{q}\right\rfloor\right). \end{aligned}$

Because \(c(p)=0\), the function \(c(q)\) vanishes whenever some prime appears to the first power only. Thus the surviving \(q\) are exactly the integers for which every prime exponent is either \(0\) or at least \(2\). These are the powerful numbers.

Step 4: Compute \(D(x)\) with the Hyperbola Method

The divisor summatory function can be written as

$D(x)=\sum_{m \le x}\tau(m)=\sum_{uv \le x}1.$

Counting lattice points under the hyperbola \(uv=x\) gives the classical identity

$D(x)=2\sum_{u=1}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\frac{x}{u}\right\rfloor-\lfloor\sqrt{x}\rfloor^2.$

This is already much faster than summing \(\tau(m)\) one value at a time. The implementation then improves it further by grouping together consecutive indices that share the same quotient \(\left\lfloor x/u\right\rfloor\), so each call runs in roughly square-root time.

Step 5: Why the Recursion Enumerates Each Powerful Number Once

Every powerful number has a unique factorization

$q=\prod_{i=1}^{r} p_i^{e_i}\qquad (e_i\ge 2).$

The recursion chooses primes in increasing order, and for each chosen prime it tries the exponents \(2,3,4,\dots\) until the product would exceed \(N\). This guarantees uniqueness: the same powerful number cannot be generated in two different branches.

The local weight also matches the definition of \(c\). If a prime already appears with exponent \(e\ge 2\), increasing that exponent by \(1\) multiplies the contribution by

$\frac{c(p^{e+1})}{c(p^e)}=\frac{2^{e-1}}{2^{e-2}}=2.$

That is why the recursive coefficient simply doubles whenever the search extends one more power of the same prime.

Worked Example: \(N=10\)

Directly,

$\begin{aligned} S(10) &= 2^{\Omega(1)}+2^{\Omega(2)}+2^{\Omega(3)}+2^{\Omega(4)}+2^{\Omega(5)} \\ &\quad +2^{\Omega(6)}+2^{\Omega(7)}+2^{\Omega(8)}+2^{\Omega(9)}+2^{\Omega(10)} \\ &= 1+2+2+4+2+4+2+8+4+4=33. \end{aligned}$

Now apply the decomposition. First,

$D(10)=\sum_{m \le 10}\tau(m)=27,\qquad D(2)=3,\qquad D(1)=1.$

The powerful numbers not exceeding \(10\) are

$1,\ 4,\ 8,\ 9,$

with weights

$c(1)=1,\qquad c(4)=1,\qquad c(8)=2,\qquad c(9)=1.$

Hence

$\begin{aligned} S(10) &= c(1)D(10)+c(4)D(2)+c(8)D(1)+c(9)D(1) \\ &= 1\cdot 27 + 1\cdot 3 + 2\cdot 1 + 1\cdot 1 \\ &= 33. \end{aligned}$

This small example shows exactly what the full algorithm does at large scale: the dense part is handled by \(D(x)\), and the sparse correction comes from powerful-number factors.

How the Code Works

The C++, Python, and Java implementations begin by generating all primes up to \(\lfloor\sqrt{N}\rfloor\). No larger prime can appear in a powerful factor \(q\le N\) with exponent at least \(2\), so this prime list is sufficient for the entire recursive search.

Next, the implementation memoizes evaluations of \(D(x)\). Each fresh query uses the hyperbola formula above, together with quotient blocks where many consecutive divisors give the same floor value. Because the recursion repeatedly asks for \(D(\lfloor N/q\rfloor)\) at overlapping arguments, this cache removes a large amount of duplicate work.

Finally, the recursive search walks through powerful numbers \(q\). Every state contributes \(D(\lfloor N/q\rfloor)\), then branches by appending a new prime with exponent \(2\) or by increasing that exponent further while the product stays within range. The accumulated sum is stored in wide integer arithmetic so that the final value is preserved exactly.

Complexity Analysis

Let \(R=\lfloor\sqrt{N}\rfloor\). Building the prime sieve costs \(O(R\log\log R)\) time and \(O(R)\) memory. A fresh computation of \(D(x)\) needs about \(O(\sqrt{x})\) quotient blocks, and memoization reuses repeated arguments across the recursion. The recursive states correspond to powerful numbers \(q\le N\), which form a sparse subset of the integers, so the search tree is far smaller than a full scan up to \(N\). In practice the overall method is sublinear in \(N\) and easily fast enough for \(N=10^{14}\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=708
2. Prime omega function: Wikipedia - Prime omega function
3. Divisor function: Wikipedia - Divisor function
4. Divisor summatory function: Wikipedia - Divisor summatory function
5. Dirichlet convolution: Wikipedia - Dirichlet convolution
6. Powerful number: Wikipedia - Powerful number