Project Euler

A Huge Binomial Coefficient

Let n = 10^18 and k = 10^9. For each triple of distinct primes (p_1, p_2, p_3) with 1000 < p_i < 5000, compute C(n, k) mod (p_1 p_2 p_3) using the Chinese Remainder Theorem. Find the sum of all suc...

Source sync Apr 19, 2026

Problem #0365

Level Level 12

Solved By 1,498

Languages C++, Python

Answer 162619462356610313

Length 281 words

modular_arithmeticnumber_theorylinear_algebra

The binomial coefficient $\displaystyle {\binom {10^{18}}{10^9}}$ is a number with more than $9$ billion ($9\times 10^9$) digits.

Let $M(n,k,m)$ denote the binomial coefficient $\displaystyle {\binom {n}{k}}$ modulo $m$.

Calculate $\displaystyle {\sum M(10^{18},10^9,p\cdot q\cdot r)}$ for $1000 < p < q < r < 5000$ and $p$,$q$,$r$ prime.

Problem 365: A Huge Binomial Coefficient

Mathematical Foundation

Theorem 1 (Lucas’ Theorem). Let $p$ be a prime and let $n, k$ be non-negative integers with base- $p$ representations $n = \sum_{i=0}^{m} n_i p^i$ and $k = \sum_{i=0}^{m} k_i p^i$ (padding with leading zeros). Then

\binom{n}{k} \equiv \prod_{i=0}^{m} \binom{n_i}{k_i} \pmod{p}

where $\binom{n_i}{k_i} = 0$ if $k_i > n_i$ .

Proof. Consider the polynomial identity in $\mathbb{F}_p[x]$ :

(1+x)^p \equiv 1 + x^p \pmod{p}

which holds because $\binom{p}{j} \equiv 0 \pmod{p}$ for $1 \le j \le p-1$ . By induction, $(1+x)^{p^i} \equiv 1 + x^{p^i} \pmod{p}$ . Therefore:

(1+x)^n = \prod_{i=0}^{m} (1+x)^{n_i p^i} = \prod_{i=0}^{m} \left((1+x^{p^i})\right)^{n_i} \pmod{p}.

Expanding the left side, the coefficient of $x^k$ is $\binom{n}{k}$ . Expanding the right side, the coefficient of $x^k = x^{\sum k_i p^i}$ is $\prod_i \binom{n_i}{k_i}$ . $\square$

Theorem 2 (Chinese Remainder Theorem). Let $m_1, m_2, \ldots, m_r$ be pairwise coprime positive integers, and let $M = \prod_{i=1}^{r} m_i$ . For any integers $a_1, \ldots, a_r$ , there exists a unique $x$ with $0 \le x < M$ such that $x \equiv a_i \pmod{m_i}$ for all $i$ . Explicitly,

x = \sum_{i=1}^{r} a_i \cdot M_i \cdot (M_i^{-1} \bmod m_i) \bmod M

where $M_i = M / m_i$ .

Proof. Existence follows from the surjectivity of the map $\mathbb{Z}/M\mathbb{Z} \to \prod_{i} \mathbb{Z}/m_i\mathbb{Z}$ , which can be verified by construction: since $\gcd(M_i, m_i) = 1$ , the inverse $M_i^{-1} \bmod m_i$ exists, and $a_i M_i (M_i^{-1} \bmod m_i) \equiv a_i \pmod{m_i}$ while $\equiv 0 \pmod{m_j}$ for $j \neq i$ . Uniqueness follows since $|{\mathbb{Z}/M\mathbb{Z}}| = \prod |{\mathbb{Z}/m_i\mathbb{Z}}|$ . $\square$

Lemma 1 (Prime Count). There are exactly 574 primes in the interval $(1000, 5000)$ . The number of unordered triples is $\binom{574}{3} = 31{,}436{,}024$ . $\square$

Lemma 2 (Base- $p$ Digits of $10^{18}$ and $10^9$ ). For a prime $p$ in $(1000, 5000)$ , the number of base- $p$ digits of $10^{18}$ is $\lfloor \log_p(10^{18}) \rfloor + 1 \le \lfloor 18 \log_{1009}(10) \rfloor + 1 \le 7$ . Computing the base- $p$ representation and evaluating the product in Lucas’ theorem takes $O(\log_p N)$ multiplications modulo $p$ , hence $O(1)$ per prime. $\square$

Editorial

Compute C(10^18, 10^9) mod p for primes p in (1000, 5000) using Lucas’ theorem. Then sum C(10^18, 10^9) mod (p1p2p3) over all triples using CRT. We compute r[p] = C(n, k) mod p for each prime p. We then iterate over p in primes. Finally, iterate over each triple, apply CRT and accumulate.

Pseudocode

Compute r[p] = C(n, k) mod p for each prime p
for p in primes
For each triple, apply CRT and accumulate

Complexity Analysis

Time: $O(\pi^3 / 6)$ for iterating over all triples, where $\pi = 574$ . This gives approximately $31.4 \times 10^6$ CRT evaluations, each $O(1)$ . Precomputing the residues takes $O(\pi \cdot \log_p N) = O(\pi)$ . Total: $O(\pi^3)$ .
Space: $O(\pi)$ for storing the per-prime residues.

Answer

\boxed{162619462356610313}

C++ project_euler/problem_365/solution.cpp

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

/*
 * Problem 365: A Huge Binomial Coefficient
 *
 * Compute C(10^18, 10^9) mod p for primes p in (1000, 5000) using Lucas' theorem.
 * Then sum C(10^18, 10^9) mod (p1*p2*p3) over all triples using CRT.
 *
 * Answer: 162619462356610313
 */

typedef long long ll;
typedef unsigned long long ull;
typedef __int128 lll;

// Sieve primes in range (lo, hi)
vector<int> sieve_primes(int lo, int hi) {
    vector<bool> is_prime(hi + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i * i <= hi; i++) {
        if (is_prime[i]) {
            for (int j = i * i; j <= hi; j += i)
                is_prime[j] = false;
        }
    }
    vector<int> primes;
    for (int i = lo + 1; i < hi; i++) {
        if (is_prime[i]) primes.push_back(i);
    }
    return primes;
}

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

// C(n, k) mod p for n, k < p (direct computation)
ll small_comb(ll n, ll k, ll p) {
    if (k > n) return 0;
    if (k == 0 || k == n) return 1;
    if (k > n - k) k = n - k;

    ll num = 1, den = 1;
    for (ll i = 0; i < k; i++) {
        num = num * ((n - i) % p) % p;
        den = den * ((i + 1) % p) % p;
    }
    return num % p * power(den, p - 2, p) % p;
}

// Lucas' theorem: C(n, k) mod p
ll lucas(ll n, ll k, ll p) {
    ll result = 1;
    while (n > 0 || k > 0) {
        ll ni = n % p, ki = k % p;
        if (ki > ni) return 0;
        result = result * small_comb(ni, ki, p) % p;
        n /= p;
        k /= p;
    }
    return result;
}

// Extended GCD
ll extgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = extgcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

// CRT for three residues
ll crt3(ll r1, ll p1, ll r2, ll p2, ll r3, ll p3) {
    lll M = (lll)p1 * p2 * p3;
    lll M1 = (lll)p2 * p3;
    lll M2 = (lll)p1 * p3;
    lll M3 = (lll)p1 * p2;

    ll x, y;
    extgcd((ll)(M1 % p1), p1, x, y);
    lll inv1 = ((lll)x % p1 + p1) % p1;

    extgcd((ll)(M2 % p2), p2, x, y);
    lll inv2 = ((lll)x % p2 + p2) % p2;

    extgcd((ll)(M3 % p3), p3, x, y);
    lll inv3 = ((lll)x % p3 + p3) % p3;

    lll result = (r1 * M1 % M * inv1 % M + r2 * M2 % M * inv2 % M + r3 * M3 % M * inv3 % M) % M;
    return (ll)result;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    ll N = 1000000000000000000LL; // 10^18
    ll K = 1000000000LL;          // 10^9

    vector<int> primes = sieve_primes(1000, 5000);
    int sz = primes.size();

    // Compute C(N, K) mod p for each prime
    vector<ll> residues(sz);
    for (int i = 0; i < sz; i++) {
        residues[i] = lucas(N, K, primes[i]);
    }

    // Sum over all triples
    lll total = 0;
    for (int i = 0; i < sz; i++) {
        for (int j = i + 1; j < sz; j++) {
            for (int k = j + 1; k < sz; k++) {
                ll val = crt3(residues[i], primes[i],
                              residues[j], primes[j],
                              residues[k], primes[k]);
                total += val;
            }
        }
    }

    cout << (ll)total << endl;
    // Expected: 162619462356610313

    return 0;
}

Python project_euler/problem_365/solution.py

"""
Problem 365: A Huge Binomial Coefficient

Compute C(10^18, 10^9) mod p for primes p in (1000, 5000) using Lucas' theorem.
Then sum C(10^18, 10^9) mod (p1*p2*p3) over all triples using CRT.

Answer: 162619462356610313
"""

from itertools import combinations

def sieve_primes(lo, hi):
    """Return list of primes in (lo, hi)."""
    is_prime = [True] * (hi + 1)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(hi**0.5) + 1):
        if is_prime[i]:
            for j in range(i * i, hi + 1, i):
                is_prime[j] = False
    return [p for p in range(lo + 1, hi) if is_prime[p]]

def small_comb(n, k, p):
    """Compute C(n, k) mod p for n, k < p."""
    if k > n:
        return 0
    if k == 0 or k == n:
        return 1
    if k > n - k:
        k = n - k
    num = 1
    den = 1
    for i in range(k):
        num = num * ((n - i) % p) % p
        den = den * ((i + 1) % p) % p
    return num * pow(den, p - 2, p) % p

def lucas(n, k, p):
    """Compute C(n, k) mod p using Lucas' theorem."""
    result = 1
    while n > 0 or k > 0:
        ni, ki = n % p, k % p
        if ki > ni:
            return 0
        result = result * small_comb(ni, ki, p) % p
        n //= p
        k //= p
    return result

def crt3(r1, p1, r2, p2, r3, p3):
    """Chinese Remainder Theorem for three congruences."""
    M = p1 * p2 * p3
    M1, M2, M3 = p2 * p3, p1 * p3, p1 * p2
    inv1 = pow(M1, p1 - 2, p1)
    inv2 = pow(M2, p2 - 2, p2)
    inv3 = pow(M3, p3 - 2, p3)
    return (r1 * M1 * inv1 + r2 * M2 * inv2 + r3 * M3 * inv3) % M

def solve():
    N = 10**18
    K = 10**9

    primes = sieve_primes(1000, 5000)

    # Compute residues
    residues = {p: lucas(N, K, p) for p in primes}

    # Sum over all triples
    total = 0
    for p1, p2, p3 in combinations(primes, 3):
        val = crt3(residues[p1], p1, residues[p2], p2, residues[p3], p3)
        total += val

    return total

if __name__ == "__main__":
    answer = solve()
assert answer == 162619462356610313
print(answer)
    # Expected: 162619462356610313