Project Euler

Divisors of Binomial Product

Let B(n) = prod_(k=0)^n C(n, k), a product of binomial coefficients. Let D(n) = sum_(d | B(n)) d, the sum of the divisors of B(n). For example, B(5) = C(5, 0) x C(5, 1) x C(5, 2) x C(5, 3) x C(5, 4...

Source sync Apr 19, 2026

Problem #0650

Level Level 10

Solved By 2,025

Languages C++, Python

Answer 538319652

Length 297 words

modular_arithmeticnumber_theorylinear_algebra

Let $B(n) = \displaystyle \prod _{k=0}^n {n \choose k}$, a product of binomial coefficients.

For example, $B(5) = {5 \choose 0} \times {5 \choose 1} \times {5 \choose 2} \times {5 \choose 3} \times {5 \choose 4} \times {5 \choose 5} = 1 \times 5 \times 10 \times 10 \times 5 \times 1 = 2500$.

Let $D(n) = \displaystyle \sum _{d|B(n)} d$, the sum of the divisors of $B(n)$.

For example, the divisors of $B(5)$ are $1$, $2$, $4$, $5$, $10$, $20$, $25$, $50$, $100$, $125$, $250$, $500$, $625$, $1250$ and $2500$,

so $D(5) = 1 + 2 + 4 + 5 + 10 + 20 + 25 + 50 + 100 + 125 + 250 + 500 + 625 + 1250 + 2500 = 5467$.

Let $S(n) = \displaystyle \sum _{k=1}^n D(k)$.

You are given $S(5) = 5736$, $S(10) = 141740594713218418$ and $S(100)$ mod $1\,000\,000\,007 = 332792866$.

Find $S(20\,000)$ mod $1\,000\,000\,007$.

Problem 650: Divisors of Binomial Product

Mathematical Analysis

Step 1: Prime Factorization of B(n)

We have $B(n) = \prod_{k=0}^{n} \binom{n}{k}$ . To find the prime factorization of $B(n)$ , we need the exponent of each prime $p$ in $B(n)$ .

The exponent of prime $p$ in $\binom{n}{k}$ is $\sum_{i=1}^{\infty} \left(\lfloor n/p^i \rfloor - \lfloor k/p^i \rfloor - \lfloor (n-k)/p^i \rfloor\right)$ .

The total exponent of $p$ in $B(n) = \prod_{k=0}^{n} \binom{n}{k}$ is:

e_p(n) = \sum_{k=0}^{n} v_p\left(\binom{n}{k}\right)

Step 2: Incremental Computation

Key identity: $\prod_{k=0}^{n} \binom{n}{k} = \prod_{k=1}^{n} \frac{k^{2k-n-1} \cdot n!}{1}$ … This simplifies with careful analysis.

A more practical approach: Track the exponent of each prime $p \le n$ in $B(n)$ . We can compute $B(n)$ incrementally from $B(n-1)$ .

Note that $\binom{n}{k} = \binom{n-1}{k} \cdot \frac{n}{n-k}$ , so:

\frac{B(n)}{B(n-1)} = \prod_{k=0}^{n} \binom{n}{k} \Big/ \prod_{k=0}^{n-1} \binom{n-1}{k} = \frac{n^n}{\prod_{k=1}^{n} k} \cdot \text{(adjustments)}

Actually, the cleanest approach uses:

B(n) = \prod_{k=0}^{n} \binom{n}{k} = \frac{n!^{n+1}}{\prod_{k=0}^{n} (k!)^2}

So the exponent of $p$ in $B(n)$ is $(n+1) \cdot v_p(n!) - 2\sum_{k=0}^{n} v_p(k!)$ .

Step 3: Sum of Divisors

If $B(n) = \prod p_i^{a_i}$ , then $D(n) = \sigma(B(n)) = \prod_i \frac{p_i^{a_i+1} - 1}{p_i - 1}$ .

We compute this modulo $10^9 + 7$ using modular arithmetic and Fermat’s little theorem for the modular inverse of $(p_i - 1)$ .

Step 4: Incremental Update

From step $n-1$ to step $n$ , we update the prime factorization. The ratio $B(n)/B(n-1)$ can be computed by noting:

\frac{B(n)}{B(n-1)} = \frac{n^n}{\prod_{k=1}^{n-1} k} = \frac{n^n}{(n-1)!}

This means going from $B(n-1)$ to $B(n)$ , we multiply by $n^n / (n-1)!$ , so we add $n \cdot v_p(n)$ and subtract $v_p((n-1)!)$ for each prime $p$ .

Editorial

We use a sieve to find smallest prime factor for all numbers up to $n = 20000$ . Finally, maintain the prime factorization of $B(k)$ incrementally.

Pseudocode

Use a sieve to find smallest prime factor for all numbers up to $n = 20000$
Maintain the prime factorization of $B(k)$ incrementally
For each $k$, compute $D(k) = \prod_p \frac{p^{e_p+1}-1}{p-1} \pmod{10^9+7}$
Sum all $D(k)$ to get $S(20000) \pmod{10^9+7}$

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: $O(n^2 / \log n)$ in the worst case for tracking all prime exponents, but effectively $O(n \cdot \pi(n))$ where $\pi(n)$ is the number of primes up to $n$ .
Space: $O(n)$ for storing prime factorizations.

Answer

\boxed{538319652}

C++ project_euler/problem_650/solution.cpp

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 1000000007;

long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

long long modinv(long long a, long long mod) {
    return power(a, mod - 2, mod);
}

int main() {
    const int N = 20000;

    // Smallest prime factor sieve
    vector<int> spf(N + 1, 0);
    vector<int> primes;
    for (int i = 2; i <= N; i++) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
            for (long long j = (long long)i * i; j <= N; j += i) {
                if (spf[j] == 0) spf[j] = i;
            }
        }
    }

    // B(n) = prod_{k=0}^{n} C(n,k)
    // The exponent of prime p in B(n):
    // v_p(B(n)) = sum_{k=0}^{n} v_p(C(n,k))
    //           = sum_{k=0}^{n} [v_p(n!) - v_p(k!) - v_p((n-k)!)]
    //           = (n+1)*v_p(n!) - 2*sum_{k=0}^{n} v_p(k!)
    //
    // Let's define cum_vp(n) = sum_{k=0}^{n} v_p(k!)
    // Then v_p(B(n)) = (n+1)*v_p(n!) - 2*cum_vp(n)
    //
    // Note: v_p(k!) = v_p((k-1)!) + v_p(k)
    // And: cum_vp(n) = cum_vp(n-1) + v_p(n!)
    //
    // Incrementally from n-1 to n:
    // v_p(B(n)) = (n+1)*v_p(n!) - 2*cum_vp(n)
    // v_p(B(n-1)) = n*v_p((n-1)!) - 2*cum_vp(n-1)
    //
    // v_p(B(n)) - v_p(B(n-1)) = (n+1)*v_p(n!) - 2*cum_vp(n) - n*v_p((n-1)!) + 2*cum_vp(n-1)
    //   = (n+1)*(v_p((n-1)!) + v_p(n)) - n*v_p((n-1)!) - 2*v_p(n!)
    //   = v_p((n-1)!) + (n+1)*v_p(n) - 2*(v_p((n-1)!) + v_p(n))
    //   = v_p((n-1)!) + (n+1)*v_p(n) - 2*v_p((n-1)!) - 2*v_p(n)
    //   = (n-1)*v_p(n) - v_p((n-1)!)
    //
    // So delta = (n-1)*v_p(n) - v_p((n-1)!)

    // Precompute v_p(k) for each k using spf
    // v_p(n!) = sum_{i>=1} floor(n/p^i)

    auto vp_factorial = [](long long n, long long p) -> long long {
        long long res = 0;
        long long pk = p;
        while (pk <= n) {
            res += n / pk;
            if (pk > n / p) break; // overflow prevention
            pk *= p;
        }
        return res;
    };

    auto vp_val = [&](int n, int p) -> int {
        int cnt = 0;
        while (n % p == 0) { n /= p; cnt++; }
        return cnt;
    };

    // For each prime, track exponent in B(n)
    // Start with B(1) = C(1,0)*C(1,1) = 1, all exponents 0
    // B(0) = C(0,0) = 1

    // We'll track exponents as a vector indexed by prime index
    int np = primes.size();
    vector<long long> expo(np, 0); // expo[i] = exponent of primes[i] in B(current_n)

    // Precompute modinv for each prime
    vector<long long> inv_pm1(np);
    for (int i = 0; i < np; i++) {
        inv_pm1[i] = modinv(primes[i] - 1, MOD);
    }

    long long S = 0;
    // B(1) = 1, D(1) = 1
    S = 1; // D(1) = 1

    // Compute for n = 2 to N
    // We need to track expo from B(1)
    // B(1): all exponents are 0

    for (int n = 2; n <= N; n++) {
        // delta for prime p: (n-1)*v_p(n) - v_p((n-1)!)
        for (int i = 0; i < np && primes[i] <= n; i++) {
            int p = primes[i];
            long long vpn = vp_val(n, p);
            if (vpn > 0 || p <= n - 1) {
                long long delta = (long long)(n - 1) * vpn - vp_factorial(n - 1, p);
                expo[i] += delta;
            }
        }

        // Compute D(n)
        long long D = 1;
        for (int i = 0; i < np && primes[i] <= n; i++) {
            if (expo[i] > 0) {
                long long num = (power(primes[i], expo[i] + 1, MOD) - 1 + MOD) % MOD;
                D = D * num % MOD * inv_pm1[i] % MOD;
            }
        }

        S = (S + D) % MOD;
    }

    cout << S << endl;
    return 0;
}

Python project_euler/problem_650/solution.py

#!/usr/bin/env python3
"""Project Euler Problem 650: Divisors of Binomial Product"""

MOD = 1000000007

def solve():
    N = 20000

    # Sieve of smallest prime factors
    spf = list(range(N + 1))
    for i in range(2, int(N**0.5) + 1):
        if spf[i] == i:  # i is prime
            for j in range(i * i, N + 1, i):
                if spf[j] == j:
                    spf[j] = i

    # Collect primes
    primes = [i for i in range(2, N + 1) if spf[i] == i]
    prime_idx = {p: i for i, p in enumerate(primes)}
    np_ = len(primes)

    def vp_val(n, p):
        cnt = 0
        while n % p == 0:
            n //= p
            cnt += 1
        return cnt

    def vp_factorial(n, p):
        res = 0
        pk = p
        while pk <= n:
            res += n // pk
            pk *= p
        return res

    # Precompute modular inverses
    inv_pm1 = [pow(p - 1, MOD - 2, MOD) if p > 1 else 0 for p in primes]

    expo = [0] * np_

    S = 1  # D(1) = 1

    for n in range(2, N + 1):
        # Update exponents: delta = (n-1)*v_p(n) - v_p((n-1)!)
        for i in range(np_):
            p = primes[i]
            if p > n:
                break
            vpn = vp_val(n, p)
            delta = (n - 1) * vpn - vp_factorial(n - 1, p)
            expo[i] += delta

        # Compute D(n)
        D = 1
        for i in range(np_):
            p = primes[i]
            if p > n:
                break
            if expo[i] > 0:
                num = (pow(p, expo[i] + 1, MOD) - 1) % MOD
                D = D * num % MOD * inv_pm1[i] % MOD

        S = (S + D) % MOD

    print(S)

solve()