Problem 365: A Huge Binomial Coefficient

Detailed mathematical approach

Problem Summary

For integers \(n\ge k\ge 0\) and a modulus \(m\ge 1\), let \(M(n,k,m)\) denote \(\binom{n}{k}\) modulo \(m\), the remainder of \(\binom{n}{k}\) on division by \(m\). We must evaluate

$$\sum_{1000\lt p\lt q\lt r\lt 5000}M(10^{18},10^{9},pqr),$$

where \(p\), \(q\) and \(r\) are prime. ^[1] The statement remarks that \(\binom{10^{18}}{10^{9}}\) itself has more than \(9\times 10^{9}\) digits; Section 1 pins the figure down to \(9\,434\,294\,477\). There are \(501\) primes strictly between \(1000\) and \(5000\), from \(1009\) to \(4999\), so the sum runs over \(\binom{501}{3}=20\,833\,250\) triples, each term below \(5000^{3}\).

The whole solution in one sentence: reduce the giant coefficient modulo each of the \(501\) primes once, using Lucas's theorem to turn it into a product of binomial coefficients of base-\(p\) digits, then glue three of those residues with the Chinese remainder theorem to get the value modulo \(pqr\), and add the \(20\,833\,250\) glued values.

Mathematical Approach

1. The number itself is out of reach

With \(n=10^{18}\) and \(k=10^{9}\), write \(\binom nk=\frac{n^{k}}{k!}\prod_{i\lt k}(1-\tfrac in)\) and apply Stirling's series \(\ln k!=k\ln k-k+\tfrac12\ln(2\pi k)+\dots\): ^[8]

$$\begin{aligned}\log_{10}\binom nk\approx{}&k\log_{10}\frac nk+k\log_{10}e\\&-\tfrac12\log_{10}(2\pi k)-\frac{k^{2}}{2n\ln 10}.\end{aligned}$$

The four terms are \(9\times 10^{9}\), \(4.343\times 10^{8}\), \(-4.899\) and \(-0.217\). The last is decisive: drop it and the estimate crosses \(9434294477\), inflating the digit count by one. Stirling's truncation error is below \(1/(12k)\lt 10^{-10}\), and a forty-digit evaluation of \(\ln n!-\ln k!-\ln(n-k)!\) confirms \(\log_{10}\binom nk=9434294476.787\ldots\); a positive integer has one more digit than the floor of its base-10 logarithm, so the coefficient has exactly \(9\,434\,294\,477\) digits. Building it is pointless rather than impossible: the remainders asked for come exactly from small integers instead.

2. One composite modulus becomes three prime ones

The moduli are products \(pqr\) of three distinct primes, which is the decisive feature. Let \(B=\binom{10^{18}}{10^{9}}\) and suppose we know

$$a=B\bmod p,\quad b=B\bmod q,\quad c=B\bmod r.$$

Claim. Exactly one \(x\) with \(0\le x\lt pqr\) satisfies \(x\equiv a\pmod p\), \(x\equiv b\pmod q\) and \(x\equiv c\pmod r\), and it is \(M(10^{18},10^{9},pqr)\).

Proof. Distinct primes are pairwise coprime, so the Chinese remainder theorem supplies such an \(x\) and makes it unique in \([0,pqr)\). ^[4] Since \(B\) satisfies the same congruences, \(B-x\) is divisible by \(p\), \(q\) and \(r\), hence by their product. So \(x\equiv B\pmod{pqr}\), and lying in \([0,pqr)\) it is the remainder of \(B\). \(\blacksquare\)

This is why the problem is tractable: the \(20\,833\,250\) composite moduli reuse only \(501\) prime residues. Had a repeated prime been allowed, say \(p^{2}q\), the argument would fail at once: Lucas's theorem below needs a prime modulus, and its prime-power analogues are far more involved. ^[2]

3. Lucas's theorem

Fix a prime \(p\) and write \(n\) and \(k\) in base \(p\), with digits \(n_t,k_t\in\{0,1,\dots,p-1\}\):

$$n=\sum_{t=0}^{s}n_tp^{t},\qquad k=\sum_{t=0}^{s}k_tp^{t}.$$

Lucas's theorem. ^[2] With the convention \(\binom ab=0\) for \(b\gt a\),

$$\binom nk\equiv\prod_{t=0}^{s}\binom{n_t}{k_t}\pmod p.$$

Proof. Work with polynomials whose coefficients are reduced modulo \(p\), written \(\mathbb F_p[X]\). For \(0\lt i\lt p\) the coefficient \(\binom pi=\frac{p}{i}\binom{p-1}{i-1}\) is divisible by \(p\), since \(p\) divides the numerator and is coprime to \(i\). Hence \((1+X)^{p}=1+X^{p}\), and raising this to the \(p\)-th power repeatedly gives \((1+X)^{p^{t}}=1+X^{p^{t}}\) for all \(t\ge 0\). Therefore

$$\begin{aligned}(1+X)^{n}&=\prod_{t=0}^{s}\bigl(1+X^{p^{t}}\bigr)^{n_t}\\&=\prod_{t=0}^{s}\sum_{j_t=0}^{n_t}\binom{n_t}{j_t}X^{j_tp^{t}}.\end{aligned}$$

Expanding the product, a term contributes to \(X^{k}\) exactly when \(\sum_t j_tp^{t}=k\) with \(0\le j_t\le n_t\lt p\). Base-\(p\) representations being unique, this forces \(j_t=k_t\) and at most one term qualifies; comparing coefficients of \(X^{k}\) gives the congruence. If some \(k_t\) exceeds \(n_t\), no term qualifies and the coefficient is \(0\). \(\blacksquare\)

4. The theorem at work on small numbers

Take \(\binom{30}{12}=86493225\). In base \(7\), \(30=(4,2)\) and \(12=(1,5)\), digits read from the most significant. The last digit of \(12\) exceeds that of \(30\), as \(5\gt 2\), so Lucas predicts \(\binom{30}{12}\equiv 0\pmod 7\); indeed \(86493225=7\cdot 12356175\). In base \(11\), \(30=(2,8)\) and \(12=(1,1)\), so

$$\binom{30}{12}\equiv\binom21\binom81=16\equiv 5\pmod{11},$$

and \(86493225=11\cdot 7863020+5\). In base \(17\), \(30=(1,13)\) and \(12=(0,12)\) give \(\binom10\binom{13}{12}=13\), matching \(86493225\bmod 17=13\). Digit-sized binomial coefficients replace an eight-digit number here, and a nine-billion-digit one in the real instance.

5. How few digits the real instance has

The interval \(1000\lt p\lt 5000\) is kind to us. Since \(p\gt 1000\) we get \(p^{6}\gt 10^{18}\), so \(10^{18}\) has at most six base-\(p\) digits; the split falls at \(10^{18/5}\approx 3981\), so \(381\) primes need six digits and the remaining \(120\), from \(3989\) upwards, need five. Likewise \(p^{3}\gt 1000^{3}=10^{9}\) and \(p^{2}\le 4999^{2}\lt 10^{9}\), so \(10^{9}\) has exactly three base-\(p\) digits throughout. The Lucas product therefore has five or six factors, of which only the lowest three can differ from \(1\); the top two or three are \(\binom{n_t}{0}=1\).

For \(p=1103\), the digits of \(10^{18}\) from the least significant are \(592,1085,846,1021,576,612\) and those of \(10^{9}\) are \(346,1055,821,0,0,0\). The three nontrivial factors are

$$\begin{aligned}\binom{592}{346}&\equiv 1045,\\\binom{1085}{1055}&\equiv 809,\\\binom{846}{821}&\equiv 442,\end{aligned}$$

all modulo \(1103\), and \(1045\cdot 809\equiv 507\), then \(507\cdot 442\equiv 185\). So \(B\equiv 185\pmod{1103}\).

6. When the residue is zero, and how often

The product vanishes precisely when some base-\(p\) digit of \(k\) exceeds its match in \(n\): if none does, every factor is a quotient of factorials of numbers below \(p\), hence nonzero, and a product of nonzero residues modulo a prime is nonzero. Kummer's theorem restates this vividly: the exponent of \(p\) in \(\binom nk\) is the number of carries when \(k\) and \(n-k\) are added in base \(p\). ^[3] For \(\binom{30}{12}\) and \(p=7\), adding \(12=(1,5)\) to \(18=(2,4)\) gives \(5+4=9\ge 7\): one carry and no second, so \(7\) divides \(\binom{30}{12}\) exactly once, and indeed \(\binom{30}{12}\bmod 49=42\ne 0\). In base \(11\) it carries nowhere, matching the residue \(5\).

In the real instance carries are common: of the \(501\) primes, \(376\) give residue \(0\) and only \(125\) a nonzero one. These zeros are a shortcut in disguise: the \(\binom{376}{3}=8\,789\,000\) triples whose residues all vanish — over two fifths of them — contribute \(0\), since \(x=0\) then satisfies all three congruences.

7. Digit binomial coefficients modulo a prime

What remains is \(\binom ab\bmod p\) for \(0\le b\le a\lt p\). An inverse of \(d\) modulo \(p\) is a residue \(v\) with \(dv\equiv 1\pmod p\); multiplying by \(v\) is division by \(d\), and modulo a prime every nonzero residue has one. ^[6] Every factor of \(b!\,(a-b)!\) lies below \(p\) and is therefore nonzero, so

$$\binom ab\equiv a!\cdot(b!)^{-1}\cdot((a-b)!)^{-1}\pmod p.$$

One pass computes \(i!\bmod p\) for all \(i\lt p\). The inverse factorials need a single inversion: by Fermat's little theorem \(x^{p-1}\equiv 1\pmod p\) for \(x\) not divisible by \(p\), so \(x^{-1}\equiv x^{p-2}\), ^[6] giving \(((p-1)!)^{-1}\) from one modular exponentiation; then \(((i-1)!)^{-1}\equiv i\,(i!)^{-1}\pmod p\) walks the table downwards. So the table costs \(O(p)\) multiplications and one exponentiation.

8. Gluing three residues: the two-step lift

The existence half of the Chinese remainder theorem can be run as a mixed-radix construction, usually named after Garner: ^[5] it lifts one modulus at a time. Start from \(x_{pq}=a+pt\), congruent to \(a\) modulo \(p\) whatever \(t\) is. Requiring \(x_{pq}\equiv b\pmod q\) gives \(pt\equiv b-a\pmod q\), hence

$$t\equiv(b-a)\,p^{-1}\pmod q.$$

Taking \(0\le t\lt q\) yields \(0\le x_{pq}\le(p-1)+p(q-1)=pq-1\). Lift again with \(x=x_{pq}+pq\,u\), which preserves both earlier congruences, and demand \(x\equiv c\pmod r\):

$$u\equiv\bigl(c-x_{pq}\bigr)(pq)^{-1}\pmod r.$$

With \(0\le u\lt r\) we get \(x\le (pq-1)+pq(r-1)\), that is \(0\le x\lt pqr\), so by Section 2 this \(x\) is exactly \(M(10^{18},10^{9},pqr)\). Each lift is a fixed handful of operations: a difference normalised modulo the new prime, a multiplication by an inverse and a multiply-add; the second also reduces \(x_{pq}\) modulo \(r\) and builds \((pq)^{-1}\) from two stored prime inverses. No modular exponentiation or extended Euclidean algorithm runs inside the triple loop.

9. One inverse table serves every triple

The lifts need \(p^{-1}\bmod q\) and \((pq)^{-1}\bmod r\). The second is not a new object, because inversion is multiplicative:

$$(pq)^{-1}\equiv p^{-1}\,q^{-1}\pmod r,$$

since \((pq)(p^{-1}q^{-1})\equiv 1\pmod r\) and inverses modulo \(r\) are unique. ^[6] One table of the \(501\cdot 500=250\,500\) pairwise inverses \(p_i^{-1}\bmod p_j\) therefore covers every triple. Modulo \(13\), for instance, \(7^{-1}\equiv 2\) and \(11^{-1}\equiv 6\) give \(2\cdot 6=12\equiv 77^{-1}\), as the rule predicts.

10. A complete miniature instance

Replace \(10^{18}\) and \(10^{9}\) by \(20\) and \(8\), so every quantity can be checked against \(\binom{20}{8}=125970\). Take \(p=7\), \(q=11\), \(r=13\). In base \(7\), \(20=(2,6)\) and \(8=(1,1)\), so \(a\equiv\binom21\binom61=12\equiv 5\pmod 7\). In base \(11\), \(20=(1,9)\) and \(8=(0,8)\), so \(b\equiv\binom10\binom98=9\pmod{11}\). In base \(13\), \(20=(1,7)\) and \(8=(0,8)\), and \(8\gt 7\) forces \(c=0\).

Now lift. With \(7^{-1}\equiv 8\pmod{11}\) we get \(t\equiv(9-5)\cdot 8=32\equiv 10\pmod{11}\) and \(x_{pq}=5+7\cdot 10=75\); indeed \(75\bmod 7=5\) and \(75\bmod 11=9\). Since \(75\equiv 10\pmod{13}\) and \(77^{-1}\equiv 12\pmod{13}\),

$$u\equiv(0-10)\cdot 12\equiv 3\cdot 12=36\equiv 10\pmod{13},$$

so \(x=75+77\cdot 10=845\), which the direct reduction confirms: \(125970=125\cdot 1001+845\). The four triples from \(\{7,11,13,17\}\) give \(845+306+663+1989=3803\), the value the packaged checkpoint compares against. Over all \(\binom63=20\) triples from \(\{5,7,11,13,17,19\}\), that is \(3\lt p\lt q\lt r\lt 20\), the same procedure gives \(12210\), which exact reduction of \(125970\) reproduces term by term.

11. How large the total can become

Every term is below its modulus, so each is at most \(4987\cdot 4993\cdot 4999-1\), just under \(1.245\times 10^{11}\), and the sum is below

$$20\,833\,250\cdot 124\,475\,554\,909\lt 2.594\times 10^{18}.$$

A sharper bound is the sum of \(pqr\) over all triples, the third elementary symmetric function of the \(501\) primes, about \(5.279\times 10^{17}\). Both are far inside the signed \(64\)-bit range \(2^{63}-1\approx 9.223\times 10^{18}\), so that type suffices for the accumulator; the implementations use wider totals still.

How the Code Works

C++ runs its checkpoints first unless --skip-checkpoints is given. All three programs then sieve the primes, reduce the coefficient modulo each, build the inverse table and accumulate over triples; the C++ file ^[9] is the reference version. primes_between(1000, 5000) (primesBetween in Java) is a sieve of Eratosthenes over \([0,5000)\) crossing out multiples of \(p\) from \(p^{2}\) onwards. ^[7] Its collecting loop excludes both endpoints, so the strict inequalities hold exactly, and it returns the \(501\) primes in ascending order, making every index triple \(i\lt j\lt k\) give \(p\lt q\lt r\).

binom_mod_prime_lucas (Java: binomModPrimeLucas) implements Sections 3, 5 and 7. It builds fact and inv_fact of length \(p\) by a forward pass, one mod_pow call for \(((p-1)!)^{-1}\) and a downward pass; mod_pow is binary exponentiation, multiplying the result in on an odd exponent, squaring the base and halving the exponent. The digit loop strips one base-\(p\) digit from nn and kk per iteration, returns \(0\) once a digit of \(k\) exceeds its match in \(n\) (Section 6), and otherwise multiplies in \(\binom{n_t}{k_t}\) from one fact and two inv_fact entries. It runs while either number is nonzero, so leading digits of \(n\) meet implicit zeros of \(k\).

sum_crt_binom_over_triples (Java: sumCrtBinomOverTriples; Python inlines these stages in solve) fills residues with one Lucas call per prime, then inv[i][j] \(=p_i^{-1}\bmod p_j\) by Fermat exponentiation, then runs three nested loops. The outer two compute the shared x_pq and pq of Section 8 once per pair, leaving the innermost loop the second lift alone: diff_c, the product inv[i][kidx] * inv[j][kidx] realising \((pq)^{-1}\) by Section 9, the multiplier u, and the accumulation of x_pq + pq * u. Differences are written (b - a + q) % q so no negative intermediate reaches the modulo operator.

Only C++ carries self-checks, in run_checkpoints. binom_exact_small computes a small binomial coefficient exactly: it takes \(k'=\min(k,n-k)\) by symmetry, lists the numerator factors \(n-k'+1,\dots,n\) and cancels each divisor \(d=2,\dots,k'\) against them with repeated std::gcd, keeping every intermediate small. One check compares Lucas with \(\binom{30}{12}\bmod p\) for \(p\in\{7,11,13,17,19,23,29\}\) — five of those residues are \(0\), so only \(p=11\) and \(p=17\) test the nonzero digit product — and the other rebuilds Section 10's tiny sum over \(\{7,11,13,17\}\) by both routes.

The languages differ mainly in their integer types. C++ and Java hold residues and inverses in \(32\)-bit integers and the lifts in \(64\)-bit ones; C++ accumulates in unsigned __int128, a compiler extension ^[10] printed digit by digit because the standard streams do not know the type. Java's total is a BigInteger built by total.add(BigInteger.valueOf(xVal)); how many objects that allocates depends on the values, since valueOf(0) returns a cached zero and adding zero returns the total unchanged, covering the \(8\,789\,000\) vanishing triples. Python's integers are unbounded. Section 11 shows a signed \(64\)-bit accumulator would suffice here, so the wider totals are a margin. Note too that inv is filled for every ordered pair \(i\ne j\) while only \(i\lt j\) entries are read, so half of the \(250\,500\) inversions are unnecessary.

Complexity and Verification

Let \(L=5000\) bound the primes and \(m=501\) count them. The sieve costs \(O(L\log\log L)\) time and \(O(L)\) memory. ^[7] The residue stage builds a factorial table of size \(p\) per prime, so \(\sum_p O(p)\) operations — that sum of primes is \(1\,472\,009\) — plus \(501\) short digit loops. The inverse table costs \(m(m-1)=250\,500\) modular exponentiations, exponents of \(10\) to \(13\) bits and so at most \(23\) multiplications each, and occupies \(m^{2}=251\,001\) machine integers, about a megabyte in C++ and Java. The triple loop runs \(\binom m3=20\,833\,250\) iterations of bounded work. Overall the time is

$$O\bigl(L\log\log L+\textstyle\sum_p p+m^{2}\log L+m^{3}\bigr),$$

whose last term dominates, and the memory is \(O(L+m^{2})\). The cubic term describes these implementations rather than a proved lower bound for the problem: their loops visit every increasing index triple, including the \(8\,789\,000\) that contribute nothing.

On one arm64 laptop, with GCC 16.1 at -O2, OpenJDK 17 and CPython 3.11, C++ took about \(0.08\) seconds and \(10\) MB, Java about \(0.46\) seconds including start-up, and Python about \(3.3\) seconds and \(23\) MB; within the Python run the residues took \(0.12\) seconds and the inverse table \(0.22\). All three produced byte-identical output, and the C++ checkpoints pass.

These further checks were made for this article and are not in the published sources. Lucas's theorem matched exact binomial coefficients for all \(0\le k\le n\le 199\) and all \(p\in\{2,3,5,7,11,13\}\). The \(501\) residues were recomputed with exact digit binomial coefficients instead of the factorial-and-Fermat tables, reproducing every value, and the Kummer carry criterion agreed with the Lucas vanishing test for all \(501\) primes. Section 10's miniature gave \(12210\) both by the two-step lift and by exact reduction, and the four-triple checkpoint \(3803\) both ways. Finally the full-size sum was recomputed with the textbook formula \(x=\bigl(\sum_{i=1}^{3}a_iN_iv_i\bigr)\bmod P\), where \(P=pqr\), \(N_i=P/p_i\) and \(v_i=N_i^{-1}\bmod p_i\) — a different construction with different inverses — and agreed with all three programs.

Footnotes and References

The references supply the background facts used above; the reduction to three prime residues, the digit counts, the vanishing statistics and the size estimates are derived here.

1. Project Euler 365 — A Huge Binomial Coefficient. The official statement notes that the binomial coefficient of 10^18 over 10^9 has more than 9 billion digits, defines M(n, k, m) as that coefficient modulo m, and asks for the sum of M(10^18, 10^9, p·q·r) over primes 1000 < p < q < r < 5000.
2. Wikipedia — Lucas's theorem. States the congruence as a product over base-p digit positions, with a zero factor where a digit of the lower index exceeds its match; gives the divisibility corollary, the proof from (1 + X)^p ≡ 1 + X^p used in Section 3, and notes that prime-power analogues are harder (Granville, 1997).
3. Wikipedia — Kummer's theorem. States that the p-adic valuation of the binomial coefficient of n over m is the number of carries when m is added to n − m in base p.
4. Wikipedia — Chinese remainder theorem. For pairwise coprime moduli it asserts that a solution exists and that exactly one integer x with 0 ≤ x < N solves the system, N being the product of the moduli — the argument reused in Section 2.
5. cp-algorithms — Chinese remainder theorem. Gives the mixed-radix form x = x_1 + x_2·m_1 + x_3·m_1·m_2 + …, each x_i inside its own modulus, and attributes those coefficients to Garner.
6. Wikipedia — Modular multiplicative inverse. An inverse exists exactly when the two numbers are coprime and is then unique; for a prime modulus it equals a raised to the power m − 2. Uniqueness validates Section 9's product rule.
7. Wikipedia — Sieve of Eratosthenes. Gives O(n log log n) time and O(n) memory for the basic array version, and the optimisation of crossing out each prime's multiples from its square onwards.
8. Wikipedia — Stirling's approximation. Gives the series ln n! ∼ n ln n − n + ½ ln(2πn) + 1/(12n) − …, whose truncation error after the ½ ln(2πn) term is bounded by the first omitted term, 1/(12n) — the bound Section 1 uses.
9. C++ — Euler365.cpp. The linked immutable revision holds binom_mod_prime_lucas, sum_crt_binom_over_triples and run_checkpoints, and is the implementation record for this article.
10. GCC manual — 128-bit integers documents unsigned __int128 as a language extension on targets with a wide enough integer mode, and Java — BigInteger documents immutable arbitrary-precision integers and the add method of the Java port.