Problem 365: A Huge Binomial Coefficient
View on Project EulerProject Euler Problem 365 Solution
EulerSolve provides an optimized solution for Project Euler Problem 365, A Huge Binomial Coefficient, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Define $$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$ For every prime \(p\) with \(1000 \lt p \lt 5000\), we need the residue \(B \bmod p\). Then, for every triple of distinct primes \(p \lt q \lt r\) in that interval, we reconstruct the unique number \(x_{pqr}\) with $$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod{p},\qquad x_{pqr}\equiv B \pmod{q},\qquad x_{pqr}\equiv B \pmod{r}.$$ The required answer is the sum of all these reconstructed values. The sieve used in the code finds exactly \(501\) primes in the interval, so the outer summation runs over $$\binom{501}{3}=20833250$$ prime triples. Mathematical Approach Step 1: Lucas's Theorem Reduces the Huge Binomial For a fixed prime \(p\), write \(N\) and \(K\) in base \(p\): $$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$ Lucas's theorem gives the congruence $$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod{p}.$$ If some digit satisfies \(K_t \gt N_t\), then \(\binom{N_t}{K_t}=0\), so the whole product is \(0\pmod p\). This is exactly why the implementation stops immediately and returns zero in that case. The interval \(1000 \lt p \lt 5000\) makes the digit loop very short. Because \(p \gt 1000\), the number \(10^{18}\) has at most six base-\(p\) digits, and \(10^9\) has at most three....
Detailed mathematical approach
Problem Summary
Define
$$N=10^{18},\qquad K=10^9,\qquad B=\binom{N}{K}.$$
For every prime \(p\) with \(1000 \lt p \lt 5000\), we need the residue \(B \bmod p\). Then, for every triple of distinct primes \(p \lt q \lt r\) in that interval, we reconstruct the unique number \(x_{pqr}\) with
$$0 \le x_{pqr} \lt pqr,\qquad x_{pqr}\equiv B \pmod{p},\qquad x_{pqr}\equiv B \pmod{q},\qquad x_{pqr}\equiv B \pmod{r}.$$
The required answer is the sum of all these reconstructed values. The sieve used in the code finds exactly \(501\) primes in the interval, so the outer summation runs over
$$\binom{501}{3}=20833250$$
prime triples.
Mathematical Approach
Step 1: Lucas's Theorem Reduces the Huge Binomial
For a fixed prime \(p\), write \(N\) and \(K\) in base \(p\):
$$N=\sum_{t=0}^{s} N_t p^t,\qquad K=\sum_{t=0}^{s} K_t p^t,\qquad 0 \le N_t,K_t \lt p.$$
Lucas's theorem gives the congruence
$$\binom{N}{K}\equiv \prod_{t=0}^{s}\binom{N_t}{K_t}\pmod{p}.$$
If some digit satisfies \(K_t \gt N_t\), then \(\binom{N_t}{K_t}=0\), so the whole product is \(0\pmod p\). This is exactly why the implementation stops immediately and returns zero in that case.
The interval \(1000 \lt p \lt 5000\) makes the digit loop very short. Because \(p \gt 1000\), the number \(10^{18}\) has at most six base-\(p\) digits, and \(10^9\) has at most three. So each residue \(B \bmod p\) is computed from only a handful of digit-level binomials.
Step 2: Each Digit-Level Binomial Is Small
For \(0 \le b \le a \lt p\), we evaluate
$$\binom{a}{b}=\frac{a!}{b!(a-b)!}\pmod{p}.$$
Because \(p\) is prime and \(a,b \lt p\), the denominator is invertible modulo \(p\). The code therefore precomputes
$$\text{fact}[i]=i!\pmod p,\qquad \text{invFact}[i]=(i!)^{-1}\pmod p,$$
and then uses
$$\binom{a}{b}\equiv \text{fact}[a]\cdot \text{invFact}[b]\cdot \text{invFact}[a-b]\pmod p.$$
The modular inverses come from Fermat's little theorem:
$$x^{p-1}\equiv 1\pmod p\quad\Longrightarrow\quad x^{-1}\equiv x^{p-2}\pmod p,$$
which is why all three implementations contain a fast modular exponentiation routine.
Worked Lucas Checkpoint
The C++ code verifies Lucas's theorem on the smaller example \(\binom{30}{12}\). Take \(p=11\). In base \(11\),
$$30=2\cdot 11+8,\qquad 12=1\cdot 11+1.$$
Lucas gives
$$\binom{30}{12}\equiv \binom{8}{1}\binom{2}{1}=8\cdot 2=16\equiv 5\pmod{11}.$$
The exact value is \(\binom{30}{12}=86493225\), and indeed \(86493225 \equiv 5 \pmod{11}\). The checkpoint in the local source repeats this comparison for several primes.
Step 3: Reconstruct One Triple by the Chinese Remainder Theorem
Suppose we already know
$$a\equiv B\pmod p,\qquad b\equiv B\pmod q,\qquad c\equiv B\pmod r,$$
with \(p,q,r\) distinct primes. Since these moduli are pairwise coprime, the Chinese Remainder Theorem guarantees a unique residue modulo \(pqr\).
The code uses a Garner-style two-stage reconstruction. First write
$$x_{pq}=a+p\,t.$$
To make this also congruent to \(b\pmod q\), we need
$$a+p\,t\equiv b\pmod q,$$
so
$$t\equiv (b-a)\,p^{-1}\pmod q.$$
Hence
$$x_{pq}=a+p\left((b-a)\,p^{-1}\bmod q\right),\qquad 0 \le x_{pq} \lt pq.$$
Now lift once more by writing
$$x=x_{pq}+pq\,u.$$
Imposing \(x\equiv c\pmod r\) yields
$$u\equiv (c-x_{pq})(pq)^{-1}\pmod r,$$
and therefore
$$x=x_{pq}+pq\left((c-x_{pq})(pq)^{-1}\bmod r\right),\qquad 0 \le x \lt pqr.$$
This final \(x\) is the value denoted \(x_{pqr}\) in the statement.
Step 4: Why the Precomputed Inverse Table Is Enough
The implementations store every pairwise inverse
$$\text{inv}[i][j]\equiv p_i^{-1}\pmod{p_j}.$$
Then the inverse of a product is obtained for free:
$$ (pq)^{-1}\equiv p^{-1}q^{-1}\pmod r.$$
So once the table of pairwise inverses has been built, each triple reconstruction uses only a few modular additions and multiplications. There is no extended Euclidean algorithm inside the cubic loop.
Worked CRT Checkpoint
The local C++ checkpoints also test the CRT stage on \(\binom{20}{8}=125970\) and the tiny prime set \(\{7,11,13,17\}\). For the triple \((7,11,13)\), the exact residues are
$$125970\equiv 5\pmod 7,\qquad 125970\equiv 9\pmod{11},\qquad 125970\equiv 0\pmod{13}.$$
First combine \(7\) and \(11\):
$$t\equiv (9-5)\cdot 7^{-1}\equiv 4\cdot 8\equiv 10\pmod{11},$$
so
$$x_{7,11}=5+7\cdot 10=75.$$
Next combine with \(13\): since \(75\equiv 10\pmod{13}\) and \(77\equiv 12\pmod{13}\),
$$u\equiv (0-10)\cdot 12^{-1}\equiv 3\cdot 12\equiv 10\pmod{13},$$
which gives
$$x_{7,11,13}=75+77\cdot 10=845.$$
This matches the direct reduction \(125970\bmod(7\cdot 11\cdot 13)=845\). Summing the four tiny triples \((7,11,13)\), \((7,11,17)\), \((7,13,17)\), and \((11,13,17)\) gives the checkpoint total \(3803\) used by the source.
Step 5: Final Summation
If \(\mathcal P\) denotes the set of primes between \(1000\) and \(5000\), the target quantity is
$$\boxed{S=\sum_{p \lt q \lt r,\; p,q,r\in\mathcal P} x_{pqr}.}$$
Running the supplied implementations yields
$$S=162619462356610313.$$
How the Code Works
The three language versions follow the same structure. They first generate the prime list with a sieve, then compute every Lucas residue \(B \bmod p\), then precompute all pairwise inverses \(p_i^{-1}\bmod p_j\), and finally iterate over all prime triples. The C++ version accumulates the sum in unsigned __int128, the Java version uses BigInteger, and the Python version relies on Python's built-in arbitrary-precision integers.
Complexity Analysis
Let \(m=501\) be the number of primes in the interval. Sieve construction up to \(5000\) is negligible, roughly \(O(5000\log\log 5000)\). The Lucas preprocessing costs
$$\sum_{p\in\mathcal P} O(p+\log_p N),$$
because each prime builds factorial and inverse-factorial tables of size \(p\), while the digit loop is tiny. Precomputing pairwise inverses costs \(O(m^2)\) time and memory. The dominant stage is the triple loop, which performs
$$\binom{m}{3}=\binom{501}{3}=20833250$$
constant-time CRT merges. Therefore the overall running time is \(O(m^3)\), and the memory usage is \(O(m^2)\) because of the inverse matrix.
Footnotes and References
- Problem page: https://projecteuler.net/problem=365
- Lucas's theorem: Wikipedia — Lucas's theorem
- Chinese remainder theorem: Wikipedia — Chinese remainder theorem
- Fermat's little theorem: Wikipedia — Fermat's little theorem
- Garner-style reconstruction: cp-algorithms — Chinese remainder theorem
Problem 365 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
i64 mod_pow(i64 base, i64 exp, i64 mod) {
i64 result = 1 % mod;
base %= mod;
while (exp > 0) {
if ((exp & 1LL) != 0LL) {
result = static_cast<i64>((static_cast<__int128>(result) * base) % mod);
}
base = static_cast<i64>((static_cast<__int128>(base) * base) % mod);
exp >>= 1LL;
}
return result;
}
std::string to_string_u128(u128 value) {
if (value == 0U) {
return "0";
}
std::string out;
while (value > 0U) {
const unsigned digit = static_cast<unsigned>(value % 10U);
out.push_back(static_cast<char>('0' + digit));
value /= 10U;
}
std::reverse(out.begin(), out.end());
return out;
}
std::vector<int> primes_between(const int lo_exclusive, const int hi_exclusive) {
std::vector<bool> is_prime(static_cast<std::size_t>(hi_exclusive), true);
if (hi_exclusive > 0) {
is_prime[0] = false;
}
if (hi_exclusive > 1) {
is_prime[1] = false;
}
for (int p = 2; p * p < hi_exclusive; ++p) {
if (!is_prime[static_cast<std::size_t>(p)]) {
continue;
}
for (int q = p * p; q < hi_exclusive; q += p) {
is_prime[static_cast<std::size_t>(q)] = false;
}
}
std::vector<int> primes;
for (int x = std::max(2, lo_exclusive + 1); x < hi_exclusive; ++x) {
if (is_prime[static_cast<std::size_t>(x)]) {
primes.push_back(x);
}
}
return primes;
}
int binom_mod_prime_lucas(const u64 n, const u64 k, const int p) {
std::vector<int> fact(static_cast<std::size_t>(p), 1);
for (int i = 1; i < p; ++i) {
fact[static_cast<std::size_t>(i)] =
static_cast<int>((static_cast<i64>(fact[static_cast<std::size_t>(i - 1)]) * i) % p);
}
std::vector<int> inv_fact(static_cast<std::size_t>(p), 1);
inv_fact[static_cast<std::size_t>(p - 1)] =
static_cast<int>(mod_pow(fact[static_cast<std::size_t>(p - 1)], p - 2, p));
for (int i = p - 1; i >= 1; --i) {
inv_fact[static_cast<std::size_t>(i - 1)] =
static_cast<int>((static_cast<i64>(inv_fact[static_cast<std::size_t>(i)]) * i) % p);
}
u64 nn = n;
u64 kk = k;
i64 result = 1;
while (nn > 0 || kk > 0) {
const int ni = static_cast<int>(nn % static_cast<u64>(p));
const int ki = static_cast<int>(kk % static_cast<u64>(p));
if (ki > ni) {
return 0;
}
i64 term = fact[static_cast<std::size_t>(ni)];
term = (term * inv_fact[static_cast<std::size_t>(ki)]) % p;
term = (term * inv_fact[static_cast<std::size_t>(ni - ki)]) % p;
result = (result * term) % p;
nn /= static_cast<u64>(p);
kk /= static_cast<u64>(p);
}
return static_cast<int>(result);
}
u128 sum_crt_binom_over_triples(const u64 n, const u64 k, const std::vector<int>& primes) {
const int m = static_cast<int>(primes.size());
std::vector<int> residues(static_cast<std::size_t>(m), 0);
for (int i = 0; i < m; ++i) {
residues[static_cast<std::size_t>(i)] = binom_mod_prime_lucas(n, k, primes[static_cast<std::size_t>(i)]);
}
std::vector<std::vector<int>> inv(static_cast<std::size_t>(m), std::vector<int>(static_cast<std::size_t>(m), 0));
for (int i = 0; i < m; ++i) {
for (int j = 0; j < m; ++j) {
if (i == j) {
continue;
}
const int mod = primes[static_cast<std::size_t>(j)];
inv[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] =
static_cast<int>(mod_pow(primes[static_cast<std::size_t>(i)] % mod, mod - 2, mod));
}
}
u128 total = 0;
for (int i = 0; i < m - 2; ++i) {
const i64 p = primes[static_cast<std::size_t>(i)];
const i64 a = residues[static_cast<std::size_t>(i)];
for (int j = i + 1; j < m - 1; ++j) {
const i64 q = primes[static_cast<std::size_t>(j)];
const i64 b = residues[static_cast<std::size_t>(j)];
const i64 diff_ab = (b - a + q) % q;
const i64 t = (diff_ab * inv[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)]) % q;
const i64 x_pq = a + p * t; // modulo p*q
const i64 pq = p * q;
for (int kidx = j + 1; kidx < m; ++kidx) {
const i64 r = primes[static_cast<std::size_t>(kidx)];
const i64 c = residues[static_cast<std::size_t>(kidx)];
const i64 diff_c = (c - (x_pq % r) + r) % r;
const i64 inv_pq_mod_r =
(static_cast<i64>(inv[static_cast<std::size_t>(i)][static_cast<std::size_t>(kidx)]) *
inv[static_cast<std::size_t>(j)][static_cast<std::size_t>(kidx)]) %
r;
const i64 u = (diff_c * inv_pq_mod_r) % r;
const i64 x = x_pq + pq * u; // modulo p*q*r
total += static_cast<u128>(x);
}
}
}
return total;
}
u128 binom_exact_small(const int n, const int k) {
const int kk = std::min(k, n - k);
std::vector<u64> num;
num.reserve(static_cast<std::size_t>(kk));
for (int x = n - kk + 1; x <= n; ++x) {
num.push_back(static_cast<u64>(x));
}
for (int d = 2; d <= kk; ++d) {
u64 rem = static_cast<u64>(d);
for (u64& v : num) {
const u64 g = std::gcd(v, rem);
if (g > 1U) {
v /= g;
rem /= g;
if (rem == 1U) {
break;
}
}
}
}
u128 result = 1;
for (const u64 v : num) {
result *= static_cast<u128>(v);
}
return result;
}
bool run_checkpoints() {
// Lucas check against small direct binomial modulo prime.
const u128 exact_30_12 = binom_exact_small(30, 12);
for (const int p : {7, 11, 13, 17, 19, 23, 29}) {
const int lucas = binom_mod_prime_lucas(30ULL, 12ULL, p);
const int direct = static_cast<int>(exact_30_12 % static_cast<u128>(p));
if (lucas != direct) {
std::cerr << "Checkpoint failed: Lucas mismatch for p=" << p << '\n';
return false;
}
}
// CRT+triple sum check on a tiny instance with exact arithmetic.
const std::vector<int> tiny_primes = {7, 11, 13, 17};
const u128 exact_20_8 = binom_exact_small(20, 8);
u128 direct_sum = 0;
for (int i = 0; i < static_cast<int>(tiny_primes.size()) - 2; ++i) {
for (int j = i + 1; j < static_cast<int>(tiny_primes.size()) - 1; ++j) {
for (int k = j + 1; k < static_cast<int>(tiny_primes.size()); ++k) {
const u64 mod = static_cast<u64>(tiny_primes[static_cast<std::size_t>(i)]) *
static_cast<u64>(tiny_primes[static_cast<std::size_t>(j)]) *
static_cast<u64>(tiny_primes[static_cast<std::size_t>(k)]);
direct_sum += (exact_20_8 % static_cast<u128>(mod));
}
}
}
const u128 crt_sum = sum_crt_binom_over_triples(20ULL, 8ULL, tiny_primes);
if (crt_sum != direct_sum) {
std::cerr << "Checkpoint failed: CRT triple sum mismatch\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
bool skip_checkpoints = false;
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
skip_checkpoints = true;
} else {
std::cerr << "Unknown argument: " << arg << '\n';
return 1;
}
}
if (!skip_checkpoints && !run_checkpoints()) {
return 2;
}
const std::vector<int> primes = primes_between(1000, 5000);
const u128 answer = sum_crt_binom_over_triples(1000000000000000000ULL, 1000000000ULL, primes);
std::cout << to_string_u128(answer) << '\n';
return 0;
}
Python
def solve():
N = 10**18
K = 10**9
def mod_pow(base, exp, mod):
result = 1 % mod
base %= mod
while exp > 0:
if exp & 1:
result = (result * base) % mod
base = (base * base) % mod
exp >>= 1
return result
def primes_between(lo_exc, hi_exc):
sieve = bytearray(b'\x01' * hi_exc)
sieve[0] = 0
if hi_exc > 1:
sieve[1] = 0
p = 2
while p * p < hi_exc:
if sieve[p]:
sieve[p*p::p] = bytearray(len(sieve[p*p::p]))
p += 1
return [x for x in range(max(2, lo_exc + 1), hi_exc) if sieve[x]]
def binom_mod_prime_lucas(n, k, p):
fact = [1] * p
for i in range(1, p):
fact[i] = (fact[i-1] * i) % p
inv_fact = [1] * p
inv_fact[p-1] = mod_pow(fact[p-1], p-2, p)
for i in range(p-1, 0, -1):
inv_fact[i-1] = (inv_fact[i] * i) % p
nn, kk = n, k
result = 1
while nn > 0 or kk > 0:
ni = nn % p
ki = kk % p
if ki > ni:
return 0
term = (fact[ni] * inv_fact[ki] % p) * inv_fact[ni - ki] % p
result = (result * term) % p
nn //= p
kk //= p
return result
primes = primes_between(1000, 5000)
m = len(primes)
# Precompute residues
residues = [binom_mod_prime_lucas(N, K, p) for p in primes]
# Precompute inverses
inv = [[0]*m for _ in range(m)]
for i in range(m):
for j in range(m):
if i != j:
mod = primes[j]
inv[i][j] = mod_pow(primes[i] % mod, mod - 2, mod)
# CRT triple sum
total = 0
for i in range(m - 2):
p = primes[i]
a = residues[i]
for j in range(i + 1, m - 1):
q = primes[j]
b = residues[j]
diff_ab = (b - a + q) % q
t = (diff_ab * inv[i][j]) % q
x_pq = a + p * t
pq = p * q
for kidx in range(j + 1, m):
r = primes[kidx]
c = residues[kidx]
diff_c = (c - (x_pq % r) + r) % r
inv_pq_mod_r = (inv[i][kidx] * inv[j][kidx]) % r
u = (diff_c * inv_pq_mod_r) % r
x = x_pq + pq * u
total += x
return str(total)
if __name__ == '__main__':
print(solve())
Java
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
public class Euler365 {
static long modPow(long base, long exp, long mod) {
long result = 1 % mod;
base %= mod;
while (exp > 0) {
if ((exp & 1) != 0) {
result = (result * base) % mod;
}
base = (base * base) % mod;
exp >>= 1;
}
return result;
}
static List<Integer> primesBetween(int loExclusive, int hiExclusive) {
boolean[] isPrime = new boolean[hiExclusive];
for (int i = 0; i < hiExclusive; i++)
isPrime[i] = true;
if (hiExclusive > 0)
isPrime[0] = false;
if (hiExclusive > 1)
isPrime[1] = false;
for (int p = 2; p * p < hiExclusive; p++) {
if (isPrime[p]) {
for (int q = p * p; q < hiExclusive; q += p) {
isPrime[q] = false;
}
}
}
List<Integer> primes = new ArrayList<>();
int start = Math.max(2, loExclusive + 1);
for (int x = start; x < hiExclusive; x++) {
if (isPrime[x])
primes.add(x);
}
return primes;
}
static int binomModPrimeLucas(long n, long k, int p) {
int[] fact = new int[p];
fact[0] = 1;
for (int i = 1; i < p; i++) {
fact[i] = (int) (((long) fact[i - 1] * i) % p);
}
int[] invFact = new int[p];
invFact[p - 1] = (int) modPow(fact[p - 1], p - 2, p);
for (int i = p - 1; i >= 1; i--) {
invFact[i - 1] = (int) (((long) invFact[i] * i) % p);
}
long nn = n;
long kk = k;
long result = 1;
while (nn > 0 || kk > 0) {
int ni = (int) (nn % p);
int ki = (int) (kk % p);
if (ki > ni)
return 0;
long term = fact[ni];
term = (term * invFact[ki]) % p;
term = (term * invFact[ni - ki]) % p;
result = (result * term) % p;
nn /= p;
kk /= p;
}
return (int) result;
}
static BigInteger sumCrtBinomOverTriples(long n, long k, List<Integer> primes) {
int m = primes.size();
int[] residues = new int[m];
for (int i = 0; i < m; i++) {
residues[i] = binomModPrimeLucas(n, k, primes.get(i));
}
int[][] inv = new int[m][m];
for (int i = 0; i < m; i++) {
for (int j = 0; j < m; j++) {
if (i == j)
continue;
int mod = primes.get(j);
inv[i][j] = (int) modPow(primes.get(i) % mod, mod - 2, mod);
}
}
BigInteger total = BigInteger.ZERO;
for (int i = 0; i < m - 2; i++) {
long p = primes.get(i);
long a = residues[i];
for (int j = i + 1; j < m - 1; j++) {
long q = primes.get(j);
long b = residues[j];
long diffAB = (b - a + q) % q;
long t = (diffAB * inv[i][j]) % q;
long xPq = a + p * t;
long pq = p * q;
for (int kidx = j + 1; kidx < m; kidx++) {
long r = primes.get(kidx);
long c = residues[kidx];
long diffC = (c - (xPq % r) + r) % r;
long invPqModR = ((long) inv[i][kidx] * inv[j][kidx]) % r;
long u = (diffC * invPqModR) % r;
long xVal = xPq + pq * u;
total = total.add(BigInteger.valueOf(xVal));
}
}
}
return total;
}
static String solve() {
List<Integer> primes = primesBetween(1000, 5000);
return sumCrtBinomOverTriples(1000000000000000000L, 1000000000L, primes).toString();
}
public static void main(String[] args) {
System.out.println(solve());
}
}