Problem 708: Twos Are All You Need

View on Project Euler

Project Euler Problem 708 Solution

EulerSolve provides an optimized solution for Project Euler Problem 708, Twos Are All You Need, with C++, Python, Java, and a step-by-step mathematical explanation.

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

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

Problem 708 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

u64 isqrt_u64(u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1ULL) <= x / (r + 1ULL)) ++r;
    while (r > x / r) --r;
    return r;
}

std::vector<int> sieve_primes(int limit) {
    if (limit < 2) return {};
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit + 1), 1U);
    is_prime[0] = is_prime[1] = 0U;

    for (int i = 2; 1LL * i * i <= limit; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) continue;
        for (int j = i * i; j <= limit; j += i) {
            is_prime[static_cast<std::size_t>(j)] = 0U;
        }
    }

    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(limit / std::max(1.0, std::log(static_cast<double>(limit)))));
    for (int i = 2; i <= limit; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) primes.push_back(i);
    }
    return primes;
}

std::string to_string_u128(u128 value) {
    if (value == 0) return "0";
    std::string out;
    while (value > 0) {
        out.push_back(static_cast<char>('0' + static_cast<int>(value % 10)));
        value /= 10;
    }
    std::reverse(out.begin(), out.end());
    return out;
}

class Solver708 {
  public:
    explicit Solver708(u64 n) : n_(n), primes_(sieve_primes(static_cast<int>(isqrt_u64(n)))) {
        alpha_cache_.reserve(1 << 20U);
    }

    u128 solve() { return dfs(1ULL, -1); }

  private:
    u64 n_;
    std::vector<int> primes_;
    std::unordered_map<u64, u64> alpha_cache_;

    u64 alpha(u64 n) {
        if (n <= 1ULL) return n;
        auto it = alpha_cache_.find(n);
        if (it != alpha_cache_.end()) return it->second;

        const u64 r = isqrt_u64(n);
        u128 s = 0;

        for (u64 i = 1; i <= r;) {
            const u64 q = n / i;
            const u64 j = std::min<u64>(r, n / q);
            s += static_cast<u128>(q) * static_cast<u128>(j - i + 1ULL);
            i = j + 1ULL;
        }

        const u64 ans = static_cast<u64>(2 * s - static_cast<u128>(r) * r);
        alpha_cache_.emplace(n, ans);
        return ans;
    }

    u128 dfs(u64 a, int i) {
        if (a > n_) return 0;

        u128 res = static_cast<u128>(alpha(n_ / a));

        for (int j = i + 1; j < static_cast<int>(primes_.size()); ++j) {
            const u64 p = static_cast<u64>(primes_[static_cast<std::size_t>(j)]);
            if (a > n_ / p / p) break;

            u64 pp = p * p;
            u128 w = 1;

            while (a <= n_ / pp) {
                res += w * dfs(a * pp, j);
                if (pp > n_ / p) break;
                pp *= p;
                w <<= 1U;
            }
        }
        return res;
    }
};

u64 brute_sum(int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] != 0) continue;
        spf[static_cast<std::size_t>(i)] = i;
        if (1LL * i * i <= n) {
            for (int j = i * i; j <= n; j += i) {
                if (spf[static_cast<std::size_t>(j)] == 0) spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }

    u64 sum = 1;
    for (int x = 2; x <= n; ++x) {
        int t = x;
        u64 term = 1;
        while (t > 1) {
            int p = spf[static_cast<std::size_t>(t)];
            while (t % p == 0) {
                t /= p;
                term <<= 1ULL;
            }
        }
        sum += term;
    }
    return sum;
}

void run_validations() {
    assert(static_cast<u64>(Solver708(100ULL).solve()) == brute_sum(100));
    assert(static_cast<u64>(Solver708(50'000ULL).solve()) == brute_sum(50'000));
    assert(static_cast<u64>(Solver708(100'000'000ULL).solve()) == 9'613'563'919ULL);
}

}  // namespace

int main(int argc, char** argv) {
    bool validate = true;
    for (int i = 1; i < argc; ++i) {
        std::string arg(argv[i]);
        if (arg == "--no-validate") validate = false;
    }

    if (validate) run_validations();

    constexpr u64 N = 100'000'000'000'000ULL;
    std::cout << to_string_u128(Solver708(N).solve()) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 100000000000000

    def isqrt(x):
        r = int(math.isqrt(x))
        while (r+1) <= x // (r+1): r += 1
        while r > x // r: r -= 1
        return r

    limit = isqrt(N)
    is_prime = bytearray(b'\x01'*(limit+1)); is_prime[0] = is_prime[1] = 0
    for i in range(2, int(limit**0.5)+1):
        if is_prime[i]:
            for j in range(i*i, limit+1, i): is_prime[j] = 0
    primes = [i for i in range(2, limit+1) if is_prime[i]]

    cache = {}
    def alpha(n):
        if n <= 1: return n
        if n in cache: return cache[n]
        r = isqrt(n); s = 0; i = 1
        while i <= r:
            q = n // i; j = min(r, n // q)
            s += q * (j - i + 1); i = j + 1
        ans = 2 * s - r * r
        cache[n] = ans; return ans

    def dfs(a, i):
        if a > N: return 0
        res = alpha(N // a)
        for j in range(i + 1, len(primes)):
            p = primes[j]
            if a > N // (p * p): break
            pp = p * p; w = 1
            while a <= N // pp:
                res += w * dfs(a * pp, j)
                if pp > N // p: break
                pp *= p; w <<= 1
        return res

    return str(dfs(1, -1))

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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.math.BigInteger;

public class Euler708 {

    static long isqrtU64(long x) {
        if (x < 0)
            return 0; // equivalent to unsigned behavior if larger than max
        long r = (long) Math.sqrt(x);
        while ((r + 1) <= x / (r + 1))
            ++r;
        while (r > 0 && r > x / r)
            --r;
        return r;
    }

    static List<Integer> sievePrimes(int limit) {
        List<Integer> primes = new ArrayList<>();
        if (limit < 2)
            return primes;

        byte[] isPrime = new byte[limit + 1];
        java.util.Arrays.fill(isPrime, (byte) 1);
        isPrime[0] = isPrime[1] = 0;

        for (int i = 2; i * i <= limit; ++i) {
            if (isPrime[i] != 0) {
                for (int j = i * i; j <= limit; j += i) {
                    isPrime[j] = 0;
                }
            }
        }

        for (int i = 2; i <= limit; ++i) {
            if (isPrime[i] != 0)
                primes.add(i);
        }
        return primes;
    }

    static class Solver708 {
        long n;
        List<Integer> primes;
        Map<Long, Long> alphaCache;

        Solver708(long n) {
            this.n = n;
            this.primes = sievePrimes((int) isqrtU64(n));
            this.alphaCache = new HashMap<>((int) Math.min(n, 1000000));
        }

        long alpha(long x) {
            if (x <= 1)
                return x;
            if (alphaCache.containsKey(x))
                return alphaCache.get(x);

            long r = isqrtU64(x);
            BigInteger s = BigInteger.ZERO;

            for (long i = 1; i <= r;) {
                long q = x / i;
                long j = Math.min(r, x / q);
                s = s.add(BigInteger.valueOf(q).multiply(BigInteger.valueOf(j - i + 1)));
                i = j + 1;
            }

            BigInteger ansBi = s.multiply(BigInteger.valueOf(2))
                    .subtract(BigInteger.valueOf(r).multiply(BigInteger.valueOf(r)));
            long ans = ansBi.longValue();
            alphaCache.put(x, ans);
            return ans;
        }

        BigInteger dfs(long a, int i) {
            if (a > n)
                return BigInteger.ZERO;

            BigInteger res = BigInteger.valueOf(alpha(n / a));

            for (int j = i + 1; j < primes.size(); ++j) {
                long p = primes.get(j);
                if (a > n / p / p)
                    break;

                long pp = p * p;
                BigInteger w = BigInteger.ONE;

                while (a <= n / pp) {
                    res = res.add(w.multiply(dfs(a * pp, j)));
                    if (pp > n / p)
                        break;
                    pp *= p;
                    w = w.shiftLeft(1);
                }
            }
            return res;
        }

        String solve() {
            return dfs(1, -1).toString();
        }
    }

    public static String solve() {
        long N = 100000000000000L;
        Solver708 solver = new Solver708(N);
        return solver.solve();
    }

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