Problem 608: Divisor Sums

View on Project Euler

Project Euler Problem 608 Solution

EulerSolve provides an optimized solution for Project Euler Problem 608, Divisor Sums, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must evaluate $$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$ for \(m=200!\) and \(n=10^{12}\), modulo \(10^9+7\). Here \(\sigma_0=\tau\) is the divisor-counting function. A direct double loop is hopeless, so the solution rewrites the divisor sum into an inclusion-exclusion over the primes of \(200!\), plus fast queries to the summatory divisor function. Mathematical Approach Write \(m=f!\) with \(f=200\), and rename the upper limit as \(N=n\). Let \(\mathbb{P}\) denote the set of prime numbers and define $$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$ For each \(p\in\mathcal{P}\), the exponent of \(p\) in \(f!\) is $$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$ The implementations exploit this prime-power structure very directly. Step 1: Factor the divisor sum prime by prime Every divisor of \(f!\) has the form $$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$ For a fixed integer \(k\), write \(b_p=v_p(k)\) for primes \(p\in\mathcal{P}\). Since \(\tau\) is multiplicative on prime powers, the sum over all divisors of \(f!\) separates into local contributions: $$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$ So the entire problem reduces to understanding the one-prime expression \(\sum_{a=0}^{e}(a+b+1)\)....

Detailed mathematical approach

Problem Summary

We must evaluate

$$D(m,n)=\sum_{d \mid m}\sum_{k=1}^{n}\sigma_0(kd),$$

for \(m=200!\) and \(n=10^{12}\), modulo \(10^9+7\). Here \(\sigma_0=\tau\) is the divisor-counting function. A direct double loop is hopeless, so the solution rewrites the divisor sum into an inclusion-exclusion over the primes of \(200!\), plus fast queries to the summatory divisor function.

Mathematical Approach

Write \(m=f!\) with \(f=200\), and rename the upper limit as \(N=n\). Let \(\mathbb{P}\) denote the set of prime numbers and define

$$\mathcal{P}=\{p\in\mathbb{P}: p\le f\}.$$

For each \(p\in\mathcal{P}\), the exponent of \(p\) in \(f!\) is

$$e_p=v_p(f!)=\sum_{j\ge 1}\left\lfloor\frac{f}{p^j}\right\rfloor.$$

The implementations exploit this prime-power structure very directly.

Step 1: Factor the divisor sum prime by prime

Every divisor of \(f!\) has the form

$$d=\prod_{p\in\mathcal{P}}p^{a_p},\qquad 0\le a_p\le e_p.$$

For a fixed integer \(k\), write \(b_p=v_p(k)\) for primes \(p\in\mathcal{P}\). Since \(\tau\) is multiplicative on prime powers, the sum over all divisors of \(f!\) separates into local contributions:

$$\sum_{d \mid f!}\tau(kd)=\left(\prod_{q\notin\mathcal{P}}(v_q(k)+1)\right)\prod_{p\in\mathcal{P}}\sum_{a=0}^{e_p}(a+b_p+1).$$

So the entire problem reduces to understanding the one-prime expression \(\sum_{a=0}^{e}(a+b+1)\).

Step 2: Rewrite the local sum with triangular numbers

Define the triangular-number function

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

Then

$$\sum_{a=0}^{e}(a+b+1)=(e+1)(b+1)+\frac{e(e+1)}{2}=T(e+1)(b+1)-T(e)b.$$

This identity is the key algebraic step. The term \(T(e+1)(b+1)\) keeps the current exponent of \(p\) inside \(k\), while the correction term \(T(e)b\) corresponds to removing one copy of \(p\) when \(p\mid k\).

Step 3: Expand the product by inclusion-exclusion

Now define

$$K=\prod_{p\in\mathcal{P}}T(e_p+1),\qquad \beta_p=\frac{T(e_p)}{T(e_p+1)}.$$

Choosing the correction term at a prime \(p\) introduces a factor \(-\beta_p\) and can only happen when \(p\) already divides \(k\). If \(Q\subseteq\mathcal{P}\) is the set of primes where we choose that correction, and

$$P_Q=\prod_{p\in Q}p,$$

then one copy of every prime in \(Q\) is removed from \(k\). Therefore, for each fixed \(k\),

$$\sum_{d \mid f!}\tau(kd)=K\sum_{Q\subseteq\mathcal{P},\ P_Q\mid k}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\tau\left(\frac{k}{P_Q}\right).$$

This is exactly the inclusion-exclusion that the implementations realize through a recursive traversal of prime subsets.

Step 4: Swap the order of summation

Now sum over \(1\le k\le N\). For a fixed subset \(Q\), the condition \(P_Q\mid k\) lets us write \(k=P_Qr\), so

$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)\sum_{r\le N/P_Q}\tau(r).$$

Introduce the summatory divisor function

$$S_\tau(x)=\sum_{r=1}^{x}\tau(r).$$

Then the whole problem becomes

$$D(f!,N)=K\sum_{Q\subseteq\mathcal{P}}(-1)^{|Q|}\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right).$$

If \(P_Q>N\), the corresponding term is zero, which explains the branch pruning in the recursive search.

Step 5: Evaluate the summatory divisor function quickly

The identity

$$S_\tau(x)=\sum_{d=1}^{x}\left\lfloor\frac{x}{d}\right\rfloor$$

counts factor pairs \((a,b)\) with \(ab\le x\). Splitting those lattice points across the hyperbola \(ab=x\) gives the standard Dirichlet-hyperbola formula

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

This is the large-\(x\) formula used by the implementation. Small values are handled by a precomputed prefix table.

Worked Example: \(m=3!=6\) and \(N=5\)

Here \(\mathcal{P}=\{2,3\}\) and \(e_2=e_3=1\). Hence

$$T(2)=3,\qquad T(1)=1,\qquad K=3\cdot 3=9,\qquad \beta_2=\beta_3=\frac{1}{3}.$$

The subset formula becomes

$$D(6,5)=9S_\tau(5)-3S_\tau(2)-3S_\tau(1)+S_\tau(0).$$

Now

$$S_\tau(5)=1+2+2+3+2=10,\qquad S_\tau(2)=3,\qquad S_\tau(1)=1,\qquad S_\tau(0)=0,$$

so

$$D(6,5)=9\cdot 10-3\cdot 3-3\cdot 1=78.$$

A direct check agrees:

$$\sum_{k=1}^{5}\tau(k)=10,\qquad \sum_{k=1}^{5}\tau(2k)=17,\qquad \sum_{k=1}^{5}\tau(3k)=19,\qquad \sum_{k=1}^{5}\tau(6k)=32,$$

and therefore

$$10+17+19+32=78.$$

This small case shows exactly how the inclusion-exclusion reproduces the brute-force sum.

How the Code Works

The C++, Python, and Java implementations first enumerate the primes up to \(200\) and compute each factorial exponent \(e_p\) with Legendre's formula. From those exponents they build the global factor \(K\) and the prime-specific ratios \(\beta_p\), always working modulo \(10^9+7\).

Next, they precompute \(S_\tau(x)\) for all \(x<10^6\) with a divisor sieve followed by prefix sums. For larger arguments they use the hyperbola formula above and store the results in a cache, so repeated requests for the same value are answered immediately.

The main recursive traversal walks through the primes in increasing order, maintains the current squarefree product \(P_Q\), the alternating sign, and the multiplicative weight \(K\prod_{p\in Q}\beta_p\). At each visited subset it adds the corresponding term

$$K\left(\prod_{p\in Q}\beta_p\right)S_\tau\left(\left\lfloor\frac{N}{P_Q}\right\rfloor\right)$$

with the appropriate sign. Because multiplying by another prime only increases \(P_Q\), the recursion stops as soon as the next product would exceed \(N\).

Complexity Analysis

Let \(B=10^6\). Building the small prefix table costs \(O(B\log B)\) time and \(O(B)\) memory, because every divisor updates all of its multiples. The prime list and factorial exponents up to \(200\) are tiny by comparison.

Let \(R\) be the number of squarefree products of primes \(\le 200\) that do not exceed \(N\). The subset traversal uses \(O(R)\) arithmetic steps. If \(X\) is the set of distinct large arguments passed to \(S_\tau\), then the uncached large-\(x\) cost is

$$O\left(\sum_{x\in X}\sqrt{x}\right),$$

because each such value is computed once with the hyperbola identity and then memoized. In practice, the pruning \(P_Q\le N\) and the cache are what make the computation feasible.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=608
  2. Divisor function: Wikipedia - Divisor function
  3. Divisor summatory function: Wikipedia - Divisor summatory function
  4. Legendre's formula: Wikipedia - Legendre's formula
  5. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
  6. Dirichlet hyperbola method: Wikipedia - Dirichlet hyperbola method

Problem 608 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>

namespace {

using i64 = long long;
using i128 = __int128_t;
using u64 = std::uint64_t;

constexpr i64 kMod = 1'000'000'007LL;
constexpr int kMainFact = 200;
constexpr u64 kMainN = 1'000'000'000'000ULL;
constexpr int kSmallLimit = 1'000'000;

constexpr int kCheckFact1 = 3;
constexpr u64 kCheckN1 = 100ULL;
constexpr i64 kCheckExpected1 = 3398LL;
constexpr int kCheckFact2 = 4;
constexpr u64 kCheckN2 = 1'000'000ULL;
constexpr i64 kCheckExpected2 = 268'882'292LL;

inline i64 mod_norm(i64 x) {
    x %= kMod;
    if (x < 0) x += kMod;
    return x;
}

inline i64 mod_mul(const i64 a, const i64 b) {
    return static_cast<i64>((static_cast<i128>(a) * b) % kMod);
}

i64 mod_pow(i64 base, i64 exp) {
    i64 result = 1 % kMod;
    base = mod_norm(base);
    while (exp > 0) {
        if (exp & 1LL) result = mod_mul(result, base);
        base = mod_mul(base, base);
        exp >>= 1LL;
    }
    return result;
}

std::vector<int> make_primes(const int n) {
    std::vector<unsigned char> mark(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    for (int p = 2; p <= n; ++p) {
        if (mark[static_cast<std::size_t>(p)] != 0) continue;
        primes.push_back(p);
        for (int q = p; q <= n; q += p) {
            mark[static_cast<std::size_t>(q)] = 1;
        }
    }
    return primes;
}

int exponent_in_factorial(int p, int n) {
    int e = 0;
    while (n > 0) {
        n /= p;
        e += n;
    }
    return e;
}

struct U64Hash {
    std::size_t operator()(u64 x) const noexcept {
        x += 0x9e3779b97f4a7c15ULL;
        x = (x ^ (x >> 30U)) * 0xbf58476d1ce4e5b9ULL;
        x = (x ^ (x >> 27U)) * 0x94d049bb133111ebULL;
        return static_cast<std::size_t>(x ^ (x >> 31U));
    }
};

std::vector<i64> small_prefix_tau;
std::unordered_map<u64, i64, U64Hash> large_prefix_tau_cache;

void build_small_prefix_tau() {
    small_prefix_tau.assign(static_cast<std::size_t>(kSmallLimit), 0);
    for (int d = 1; d < kSmallLimit; ++d) {
        for (int m = d; m < kSmallLimit; m += d) {
            ++small_prefix_tau[static_cast<std::size_t>(m)];
        }
    }
    for (int i = 1; i < kSmallLimit; ++i) {
        i64 v = small_prefix_tau[static_cast<std::size_t>(i)] +
                small_prefix_tau[static_cast<std::size_t>(i - 1)];
        if (v >= kMod) v -= kMod;
        small_prefix_tau[static_cast<std::size_t>(i)] = v;
    }
}

i64 prefix_tau(const u64 n) {
    if (n < static_cast<u64>(kSmallLimit)) {
        return small_prefix_tau[static_cast<std::size_t>(n)];
    }
    const auto it = large_prefix_tau_cache.find(n);
    if (it != large_prefix_tau_cache.end()) {
        return it->second;
    }

    i128 sum = 0;
    u64 x = 1;
    while (x * x <= n) {
        sum += static_cast<i128>(n / x);
        ++x;
    }
    const i128 val = 2 * sum - static_cast<i128>(x - 1) * static_cast<i128>(x - 1);
    const i64 out = mod_norm(static_cast<i64>(val % kMod));
    large_prefix_tau_cache.emplace(n, out);
    return out;
}

struct Solver608 {
    std::vector<int> primes;
    std::vector<i64> ds;
    u64 N = 0;
    i64 ans = 0;

    static i64 tri(const i64 x) { return x * (x + 1LL) / 2LL; }

    void dfs(const u64 acc, const int start_idx, const i64 d, const int sign) {
        const i64 term = mod_mul(d, prefix_tau(N / acc));
        ans = (sign > 0) ? mod_norm(ans + term) : mod_norm(ans - term);

        for (int j = start_idx; j < static_cast<int>(primes.size()); ++j) {
            const u64 p = static_cast<u64>(primes[static_cast<std::size_t>(j)]);
            if (acc > N / p) break;
            dfs(acc * p, j + 1, mod_mul(d, ds[static_cast<std::size_t>(j)]), -sign);
        }
    }

    i64 compute(const int fact_n, const u64 n) {
        N = n;
        ans = 0;
        primes = make_primes(fact_n);
        ds.assign(primes.size(), 0);

        i64 k = 1;
        for (std::size_t i = 0; i < primes.size(); ++i) {
            const int e = exponent_in_factorial(primes[i], fact_n);
            const i64 t_e = tri(e) % kMod;
            const i64 t_ep1 = tri(e + 1) % kMod;
            ds[i] = mod_mul(t_e, mod_pow(t_ep1, kMod - 2));
            k = mod_mul(k, t_ep1);
        }

        dfs(1ULL, 0, k, +1);
        return ans;
    }
};

}  // namespace

int main() {
    build_small_prefix_tau();
    large_prefix_tau_cache.reserve(1 << 16);

    Solver608 solver;

    assert(solver.compute(kCheckFact1, kCheckN1) == mod_norm(kCheckExpected1));
    assert(solver.compute(kCheckFact2, kCheckN2) == mod_norm(kCheckExpected2));

    std::cout << solver.compute(kMainFact, kMainN) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1000000007
    FACT = 200; N = 1000000000000; SL = 1000000

    def mod_norm(x):
        x %= MOD
        return x + MOD if x < 0 else x
    def mod_mul(a, b): return a * b % MOD
    def mod_pow(base, exp):
        r = 1; base = mod_norm(base)
        while exp > 0:
            if exp & 1: r = mod_mul(r, base)
            base = mod_mul(base, base); exp >>= 1
        return r

    # Sieve smallest primes up to FACT
    primes = []
    sieve = bytearray(b'\x01')*(FACT+1); sieve[0] = sieve[1] = 0
    for p in range(2, FACT+1):
        if sieve[p]:
            primes.append(p)
            for q in range(p, FACT+1, p): sieve[q] = 0

    def exp_in_fact(p, n):
        e = 0; m = n
        while m > 0: m //= p; e += m
        return e

    # Build small prefix tau
    spt = [0]*SL
    for d in range(1, SL):
        for m in range(d, SL, d): spt[m] += 1
    for i in range(1, SL): spt[i] = (spt[i] + spt[i-1]) % MOD

    cache = {}
    def prefix_tau(n):
        if n < SL: return spt[n]
        if n in cache: return cache[n]
        s = 0; x = 1
        while x*x <= n: s += n//x; x += 1
        val = mod_norm((2*s - (x-1)*(x-1)) % MOD)
        cache[n] = val; return val

    def tri(x): return x*(x+1)//2

    ds = [0]*len(primes); k = 1
    for i in range(len(primes)):
        e = exp_in_fact(primes[i], FACT)
        te = tri(e) % MOD; tep1 = tri(e+1) % MOD
        ds[i] = mod_mul(te, mod_pow(tep1, MOD-2))
        k = mod_mul(k, tep1)

    ans = [0]
    def dfs(acc, start_idx, d, sign):
        term = mod_mul(d, prefix_tau(N // acc))
        ans[0] = mod_norm(ans[0] + term * sign)
        for j in range(start_idx, len(primes)):
            p = primes[j]
            if acc > N // p: break
            dfs(acc*p, j+1, mod_mul(d, ds[j]), -sign)

    dfs(1, 0, k, 1)
    return str(ans[0])

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

Java

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

public class Euler608 {
    static final long kMod = 1000000007L;
    static final int kSmallLimit = 1000000;

    static int[] smallPrefixTau = new int[kSmallLimit];
    static Map<Long, Long> largePrefixTauCache = new HashMap<>();

    static void buildSmallPrefixTau() {
        for (int d = 1; d < kSmallLimit; d++) {
            for (int m = d; m < kSmallLimit; m += d) {
                smallPrefixTau[m]++;
            }
        }
        for (int i = 1; i < kSmallLimit; i++) {
            long v = smallPrefixTau[i] + smallPrefixTau[i - 1];
            smallPrefixTau[i] = (int) (v % kMod);
        }
    }

    static long prefixTau(long n) {
        if (n < kSmallLimit) {
            return smallPrefixTau[(int) n];
        }
        Long cached = largePrefixTauCache.get(n);
        if (cached != null)
            return cached;

        long sumVal = 0;
        long x = 1;
        while (x * x <= n) {
            sumVal += n / x;
            x++;
        }

        long xMinus1 = x - 1;
        long val = 2 * sumVal - xMinus1 * xMinus1;
        long out = val % kMod;
        if (out < 0)
            out += kMod;
        largePrefixTauCache.put(n, out);
        return out;
    }

    static List<Integer> makePrimes(int n) {
        byte[] mark = new byte[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int p = 2; p <= n; p++) {
            if (mark[p] == 0) {
                primes.add(p);
                for (int q = p; q <= n; q += p) {
                    mark[q] = 1;
                }
            }
        }
        return primes;
    }

    static int exponentInFactorial(int p, int n) {
        int e = 0;
        while (n > 0) {
            n /= p;
            e += n;
        }
        return e;
    }

    static long modPow(long base, long exp) {
        long result = 1;
        base %= kMod;
        if (base < 0)
            base += kMod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                result = (result * base) % kMod;
            base = (base * base) % kMod;
            exp >>= 1;
        }
        return result;
    }

    static long tri(long x) {
        return x * (x + 1) / 2;
    }

    static long ans = 0;
    static long N = 1000000000000L;
    static List<Integer> primes;
    static long[] ds;

    static void dfs(long acc, int startIdx, long d, int sign) {
        long term = (d * prefixTau(N / acc)) % kMod;
        if (sign > 0) {
            ans = (ans + term) % kMod;
        } else {
            ans = (ans - term + kMod) % kMod;
        }

        for (int j = startIdx; j < primes.size(); j++) {
            long p = primes.get(j);
            if (acc > N / p)
                break;
            dfs(acc * p, j + 1, (d * ds[j]) % kMod, -sign);
        }
    }

    public static String solve() {
        int factN = 200;
        buildSmallPrefixTau();

        primes = makePrimes(factN);
        ds = new long[primes.size()];

        long k = 1;
        for (int i = 0; i < primes.size(); i++) {
            int e = exponentInFactorial(primes.get(i), factN);
            long tE = tri(e) % kMod;
            long tEp1 = tri(e + 1) % kMod;
            ds[i] = (tE * modPow(tEp1, kMod - 2)) % kMod;
            k = (k * tEp1) % kMod;
        }

        ans = 0;
        dfs(1, 0, k, 1);

        return Long.toString(ans);
    }

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