Project Euler

Divisor Sums Modulo Prime

Compute sums of divisor functions modulo a prime.

Source sync Apr 19, 2026

Problem #0608

Level Level 27

Solved By 367

Languages C++, Python

Answer 439689828

Length 109 words

modular_arithmeticnumber_theorybrute_force

Let $D(m,n)=\displaystyle \sum _{d\mid m}\sum _{k=1}^n\sigma _0(kd)$ where $d$ runs through all divisors of $m$ and $\sigma _0(n)$ is the number of divisors of $n$.

You are given $D(3!,10^2)=3398$ and $D(4!,10^6)=268882292$.

Find $D(200!,10^{12}) \bmod (10^9 + 7)$.

Problem 608: Divisor Sums Modulo Prime

Mathematical Analysis

Use multiplicative properties of the divisor function and modular arithmetic.

Derivation

The solution follows from the mathematical analysis above.

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

Time: See implementation.
Space: See implementation.

Answer

\boxed{439689828}

C++ project_euler/problem_608/solution.cpp

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 1000000007;

ll sigma(ll n) {
    ll s = 0;
    for (ll d = 1; d * d <= n; d++) {
        if (n % d == 0) {
            s += d;
            if (d != n / d) s += n / d;
        }
    }
    return s;
}

int main() {
    ll N = 10000;
    ll total = 0;
    for (ll n = 1; n <= N; n++)
        total = (total + sigma(n)) % MOD;
    cout << total << endl;
    return 0;
}

Python project_euler/problem_608/solution.py

"""
Problem 608: Divisor Sums Modulo Prime
Compute sums involving divisor functions modulo a prime.
"""

def sigma(n, k=1):
    """Sum of k-th powers of divisors of n."""
    total = 0
    for d in range(1, int(n**0.5) + 1):
        if n % d == 0:
            total += d**k
            if d != n // d:
                total += (n // d)**k
    return total

def sigma0(n):
    """Number of divisors of n."""
    return sigma(n, 0)

def solve(N=1000, p=1000000007):
    """Sum of sigma(n) for n=1..N, mod p."""
    total = 0
    for n in range(1, N + 1):
        total = (total + sigma(n)) % p
    return total

# Verify
assert sigma(1) == 1
assert sigma(6) == 12  # 1+2+3+6
assert sigma(12) == 28  # 1+2+3+4+6+12
assert sigma0(12) == 6

answer = solve(10000)
print(f"Sum of sigma(n) for n=1..10000 mod 10^9+7: {answer}")