Problem 231: Prime Factorisation of Binomial Coefficients

Project Euler Problem 231 Solution

EulerSolve provides an optimized solution for Project Euler Problem 231, Prime Factorisation of Binomial Coefficients, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The target quantity is the sum of the prime factors of the binomial coefficient $\binom{20{,}000{,}000}{15{,}000{,}000}$, where each prime is counted with multiplicity. If a positive integer $N$ has prime factorisation $$N=\prod_{p \text{ prime}} p^{e_p},$$ then the required weighted sum is $$S(N)=\sum_{p \text{ prime}} p\,e_p.$$ So for this problem we must compute $$S\!\left(\binom{20{,}000{,}000}{15{,}000{,}000}\right)=\sum_{p \le 20{,}000{,}000} p\,v_p\!\left(\binom{20{,}000{,}000}{15{,}000{,}000}\right),$$ where $v_p(x)$ is the exponent of the prime $p$ in $x$. Directly expanding or factoring the binomial coefficient is impossible in practice, so the solution works entirely with prime exponents. Mathematical Approach The whole method rests on one idea: instead of building the enormous integer $\binom{n}{k}$, determine how many times each prime occurs in it and add the weighted contributions $p \cdot v_p$. Prime Exponents Inside a Binomial Coefficient Write $$\binom{n}{k}=\frac{n!}{k!(n-k)!}.$$ Prime exponents are additive across multiplication and subtractive across division, so for every prime $p$, $$v_p\!\left(\binom{n}{k}\right)=v_p(n!) - v_p(k!) - v_p((n-k)!).$$ This is the exact cancellation pattern used by the implementations. It turns the problem into three factorial-valuation queries per prime....

Detailed mathematical approach

Problem Summary

The statement writes the binomial coefficient $\binom{10}{3}=120$ as the product of primes $2\cdot2\cdot2\cdot3\cdot5$ and adds the terms of this prime factorisation: $2+2+2+3+5=14$. Each prime is counted as many times as it occurs in the factorisation. We must find the same sum for $\binom{20\,000\,000}{15\,000\,000}$. ^[1] The sum of the prime factors of $m$ with repetition is known as $\operatorname{sopfr}(m)$; ^[3] we write $f(m)$, so $f\bigl(\binom{10}{3}\bigr)=14$ is given and $f\bigl(\binom{20\,000\,000}{15\,000\,000}\bigr)$ is asked for.

The whole solution in one sentence: the binomial coefficient has almost five million digits, so instead of building it we use Legendre's formula to count, for every prime $p$ up to $20\,000\,000$, how often $p$ divides the numerator $20\,000\,000!$ and the two factorials $15\,000\,000!$ and $5\,000\,000!$ of the denominator, and add up $p$ times the difference.

Mathematical Approach

For the real problem $n=20\,000\,000$, $k=15\,000\,000$ and $r=n-k=5\,000\,000$, the names used in the programs; the small examples use the same letters for smaller values. For a prime $p$ and a positive integer $m$, $v_p(m)$ is the exponent of $p$ in $m$, the number of times $p$ divides $m$; for instance $v_2(120)=3$, $v_3(120)=1$ and $v_7(120)=0$.

1. Adding up prime factors

With this notation $f(m)=\sum_p p\,v_p(m)$, where only finitely many terms are not $0$. The empty product $1$ has $f(1)=0$; further $f(8)=2+2+2=6$, $f(12)=2+2+3=7$, and $f(p)=p$ for every prime $p$. By the fundamental theorem of arithmetic every integer greater than $1$ is a product of primes in exactly one way, up to the order of the factors, so the factorisation of $ab$ consists of those of $a$ and $b$ side by side. Hence $v_p(ab)=v_p(a)+v_p(b)$ for every prime, ^[2] and $f(ab)=f(a)+f(b)$: $f$ is totally additive. ^[3] Two consequences are used below: $f(a/b)=f(a)-f(b)$ whenever $b$ divides $a$, and $f(m!)=f(1)+f(2)+\dots+f(m)$, for example

$$\begin{aligned}f(10!)&=0+2+3+4+5\\&\quad+5+7+6+6+7=45.\end{aligned}$$

2. From one huge number to three factorials

The binomial coefficient is $\binom nk=\frac{n!}{k!\,r!}$. ^[4] For our $n$ and $k$ it has $4\,884\,378$ decimal digits, and dividing it by the primes again and again would be very slow. We do not even have to build it: since $n!=\binom nk\cdot k!\cdot r!$, the exponents subtract, and the exponent $e_p=v_p\bigl(\binom nk\bigr)$ of a prime $p$ in the binomial coefficient is

$$e_p=v_p(n!)-v_p(k!)-v_p(r!).$$

The prime factors of $n!=1\cdot2\cdots n$ are those of the numbers $2,\dots,n$, so only the primes $p\le n$ can have $e_p\gt0$, and

$$f\Bigl(\binom nk\Bigr)=\sum_{p\le n}p\,e_p.$$

Equivalently, $f\bigl(\binom nk\bigr)=f(n!)-f(k!)-f(r!)$; for the example of the statement this is $f(10!)-f(3!)-f(7!)$${}=45-5-26=14$.

3. Legendre's formula

Write $\lfloor x\rfloor$ for the largest integer not exceeding $x$; for example $\lfloor10/3\rfloor=3$. The exponent of a prime in a factorial needs no factorisation at all: ^[5]

$$v_p(m!)=\sum_{j\ge1}\Bigl\lfloor\frac{m}{p^j}\Bigr\rfloor.$$

All terms with $p^j\gt m$ are $0$, so the sum is finite. To prove the formula, write $v_p(m!)=v_p(1)+v_p(2)$${}+\dots+v_p(m)$ by Section 1. A number $i\le m$ contributes one unit for each $j\ge1$ such that $p^j$ divides $i$; $8$ contributes three units to $v_2(10!)$, for $2$, $4$ and $8$. Collecting the units by $j$ instead of by $i$, the power $p^j$ collects one unit from each of its $\lfloor m/p^j\rfloor$ multiples $p^j, 2p^j, \dots$ up to $m$. Up to $10$ there are five multiples of $2$, two of $4$ and one of $8$, so $v_2(10!)=5+2+1=8$, and indeed $10!=3628800=2^8\cdot3^4\cdot5^2\cdot7$.

Since $\bigl\lfloor\lfloor m/p^j\rfloor/p\bigr\rfloor=\lfloor m/p^{j+1}\rfloor$, ^[6] each term is the integer quotient of the previous one by $p$: for $m=10$ and $p=2$ the quotients are $5, 2, 1, 0$. The programs compute the exponent this way, with one integer division per term and without ever forming the power $p^j$.

4. Worked example: $\binom{10}{3}$

Here $n=10$, $k=3$ and $r=7$, and the primes up to $10$ are $2$, $3$, $5$ and $7$. Legendre's formula gives:

$p=2$: $v_2(10!)=5+2+1=8$, $v_2(3!)=1$ and $v_2(7!)=3+1=4$, so $e_2=8-1-4=3$;
$p=3$: $v_3(10!)=3+1=4$, $v_3(3!)=1$ and $v_3(7!)=2$, so $e_3=4-1-2=1$;
$p=5$: $v_5(10!)=2$, $v_5(3!)=0$ and $v_5(7!)=1$, so $e_5=2-0-1=1$;
$p=7$: $v_7(10!)=1$, $v_7(3!)=0$ and $v_7(7!)=1$, so $e_7=1-0-1=0$.

So $\binom{10}{3}=2^3\cdot3\cdot5=120$, and the sum is $2\cdot3+3\cdot1+5\cdot1+7\cdot0=14$, the value of the statement. The prime $7$ divides the numerator $10!$ once and the denominator $3!\cdot7!$ once, so it cancels: a prime up to $n$ need not divide the binomial coefficient.

5. Each power of $p$ adds $0$ or $1$

By Legendre's formula, $e_p$ has one term for each power $q=p^j$:

$$\begin{aligned}e_p&=\sum_{j\ge1}c_{p^j},\\c_q&=\Bigl\lfloor\frac nq\Bigr\rfloor-\Bigl\lfloor\frac kq\Bigr\rfloor-\Bigl\lfloor\frac rq\Bigr\rfloor.\end{aligned}$$

Claim. Every term $c_q$ is $0$ or $1$, and $c_q=1$ exactly when the remainders of $k$ and $r$ after division by $q$ add up to at least $q$.

Proof. Divide with remainder: $k=aq+s$ and $r=bq+t$ with $0\le s,t\lt q$. Then $n=k+r=(a+b)q+(s+t)$ with $0\le s+t\le 2q-2$. If $s+t\lt q$, then $\lfloor n/q\rfloor=a+b$ and $c_q=0$. If $s+t\ge q$, then $n=(a+b+1)q+(s+t-q)$ with $0\le s+t-q\lt q$, so $\lfloor n/q\rfloor=a+b+1$ and $c_q=1$. ^[6] $\square$

So $e_p\ge0$, and the subtraction in the programs never goes below zero. Since $c_q=0$ for $q\gt n$, $e_p$ is at most the number of powers $p^j\le n$, that is, $p^{e_p}\le n$; a prime above $\sqrt n$ divides $\binom nk$ at most once. Finally, the remainder of $k$ after division by $p^j$ is the number formed by the last $j$ digits of $k$ in base $p$, so $c_{p^j}=1$ exactly when adding $k$ and $r$ in base $p$ carries out of the $j$-th digit from the right. Thus $e_p$ is the number of carries, which is Kummer's theorem. ^[7] For example $3=11_2$ and $7=111_2$ add up to $1010_2$ with three carries, and $e_2=3$ in Section 4.

The scaled-down problem $\binom{20}{15}=15504=2^4\cdot3\cdot17\cdot19$, with $r=5$, shows this at work. For $p=17$ and $p=19$ the remainders $15$ and $5$ add up to $20\ge p$, while for $p=11$ and $p=13$ the sums $4+5$ and $2+5$ stay below $p$. For $p=2$ all of $q=2,4,8,16$ give a carry, since the sums $1+1$, $3+1$, $7+5$ and $15+5$ reach $q$, so $e_2=4$; for $p=3$ only $q=9$ does, as $6+5\ge9$, and $5$ and $7$ give none. The sum is $2\cdot4+3+17+19=47$. In the real problem the same test settles large primes in blocks: each of the $299\,903$ primes with $k\lt p\le n$ appears once, as the remainders $k$ and $r$ add up to $n\ge p$, and the $306\,125$ primes with $n/2\lt p\le k$ cancel, as $k-p$ and $r$ add up to $n-p\lt p$.

6. Listing the primes: the sieve of Eratosthenes

Section 2 needs every prime up to $n$. The sieve of Eratosthenes ^[8] marks $2,\dots,n$ and goes through $i=2,3,\dots$ while $i^2\le n$; whenever $i$ is still marked, it crosses out $i^2, i^2+i, i^2+2i, \dots$ up to $n$. A prime is never crossed out, because every crossed-out number is a product of two integers that are at least $2$. A composite $m\le n$ with smallest prime factor $q$ has $m/q\ge q$, because all prime factors of $m/q$ are at least $q$, so $q^2\le m\le n$: the loop reaches $i=q$ while $q$ is still marked and crosses out $m$. For $n=20\,000\,000$ the $607$ primes up to $4472$ make $46\,588\,036$ crossing-out steps, about $2.3$ per number, and $1\,270\,607$ primes remain.

7. The size of the sum

For $a,b\ge2$ we have $ab-(a+b)$${}=(a-1)(b-1)-1\ge0$, so a product of primes is at least their sum, and $f(m)\le m$ for every $m\ge1$, with equality only for $m=4$ and the primes. ^[3] Since $f(k!)$ and $f(r!)$ are not negative, Sections 1 and 2 give

$$\begin{aligned}f\Bigl(\binom nk\Bigr)&\le f(n!)=\sum_{m=1}^{n}f(m)\\&\le\frac{n(n+1)}2,\end{aligned}$$

about $2\cdot10^{14}$ for our $n$, far below $2^{63}\approx9.2\cdot10^{18}$. The contributions $p\,e_p$ are not negative, so every partial sum obeys the same bound, and 64-bit integers are wide enough.

8. Other ways to the same sum

Section 1 also gives a method without Legendre's formula. By the symmetry $\binom nk=\binom nr$ and the multiplicative formula, ^[4] $\binom nk$ is the product $(k+1)(k+2)\cdots n$ divided by $r!$, so

$$f\Bigl(\binom nk\Bigr)=\sum_{m=k+1}^{n}f(m)-\sum_{m=2}^{r}f(m).$$

For $\binom{20}{15}$ this is $(8+17+8+19+9)$${}-(2+3+4+5)$${}=61-14=47$. A sieve that stores the smallest prime factor $q$ of every $m\le n$ yields all values through $f(m)=q+f(m/q)$, at the cost of a table entry for every number; we use it below as an independent check. The digit-sum form of Legendre's formula ^[5] gives a third way: with $s_p(m)$ the sum of the base-$p$ digits of $m$, ^[7]

$$e_p=\frac{s_p(k)+s_p(r)-s_p(n)}{p-1}.$$

How the Code Works

The C++ program ^[9] keeps $n$, $k$ and a checkpoint switch in Options. parse_arguments accepts --n=, --k= and --skip-checkpoints, and an unknown argument or values outside $1\le n$, $0\le k\le n$ give exit code $1$. Its digit reader parse_int_after_prefix has no overflow check, so an oversized value is not rejected; a test build read --k=4294967299 as $3$. Unless it is skipped, run_checkpoints requires solve(10, 3) to return the $14$ of the statement, or the program stops with exit code $2$; then main prints the result of solve.

primes_up_to is the sieve of Section 6 on a std::vector<std::uint8_t>, one byte per number, and collects the primes into a std::vector<int>. exponent_in_factorial divides x by p until it reaches $0$ and adds up the quotients in a 64-bit counter, as in Section 3. solve sets r = n - k, forms $e_p$ for every prime from three such calls and adds $p\,e_p$ to the unsigned 64-bit sum. A negative difference would wrap around to a huge number, which Section 5 rules out, and Section 7 shows that the sum fits.

Python and Java have the same three functions, in Java primesUpTo, exponentInFactorial and solve, with the parameters as default arguments or constants and without options or checkpoint. Python sieves a bytearray with the outer loop range(2, int(n**0.5) + 1), whose floating-point bound is $4473$ here, so $i$ runs up to $4472$ as in C++. Its integers have unlimited precision, ^[9] and solve returns a string. Java sieves a byte[], counts the primes to allocate an int[] of the exact size, fills it in a second pass and adds (long) p * exp in a long.

Complexity and Verification

The sieve takes $O(n\log\log n)$ steps, ^[8] here $46\,588\,036$ crossing-out steps, followed by one pass over all numbers to collect the primes, two in Java. The flags take $n+1$ bytes, about 20 MB; the $1\,270\,607$ primes take about 5 MB as 32-bit integers in C++ and Java and about 46 MB as a Python list of integer objects. exponent_in_factorial makes $\lfloor\log_p m\rfloor+1$ divisions for $m\ge1$, so the three calls need at most $3(\lfloor\log_2 n\rfloor+1)$ divisions, here $75$, for each of the fewer than $\sqrt n$ primes up to $\sqrt n$, and at most $6$ for each larger prime: $O(\sqrt n\log n+\pi(n))=O(n)$ divisions in all, where $\pi(n)$ counts the primes up to $n$, and $6\,403\,458$ here. The sieve dominates: the method takes $O(n\log\log n)$ time and $O(n)$ memory. On an Apple M4 Max laptop, measured with /usr/bin/time, a full run takes about $0.04$ s in C++ (-O2, checkpoint included), $0.09$ s in Java and $1.25$ s in Python, with a peak memory of about 34, 69 and 83 MB.

The published code checks only $f\bigl(\binom{10}{3}\bigr)=14$, in C++. For this article the C++ program, run with --n= and --k=, agreed on all $20\,300$ pairs with $1\le n\le200$ and $0\le k\le n$ with a brute force that builds $\binom nk$ as an exact big integer and divides out every prime. Both give $47$, $2593$, $167155$ and $12953969$ for the scaled-down problems $\binom{20}{15}$, $\binom{200}{150}$, $\binom{2000}{1500}$ and $\binom{20000}{15000}$. The method of Section 8, written separately in C++, agrees with the program on $100$ random pairs with $n\le10^6$. At the real parameters the three published programs, this separate program and a Python count of the carries of Section 5 give the same number. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts; the bound of Section 7 and the operation counts are derived in this article.

Project Euler 231 — Prime Factorisation of Binomial Coefficients. The official statement factorises C(10, 3) = 120 = 2 × 2 × 2 × 3 × 5, adds the terms to 14, and asks for the sum of the terms in the prime factorisation of C(20 000 000, 15 000 000).
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a prime or a product of primes, unique up to the order of the factors; in a product, the exponents of each prime add.
OEIS A001414 — Integer log of n: sum of primes dividing n (with repetition). The sequence 0, 2, 3, 4, 5, 5, 7, 6, 6, 7, 11, 7, …, also called sopfr(n), is totally additive with a(pᵉ) = p·e, and a(n) ≤ n, with equality exactly when n is 4 or a prime.
Wikipedia — Binomial coefficient. C(n, k) = n!/(k!(n − k)!), the multiplicative formula n(n − 1)⋯(n − k + 1)/k! and the symmetry C(n, k) = C(n, n − k).
Wikipedia — Legendre's formula. The exponent of a prime p in n! is the sum of ⌊n/pⁱ⌋ over i ≥ 1 and equals (n − sₚ(n))/(p − 1), where sₚ(n) is the sum of the base-p digits of n.
Wikipedia — Floor and ceiling functions. Nested divisions satisfy ⌊⌊x/m⌋/n⌋ = ⌊x/(mn)⌋ for a positive integer n, and ⌊x⌋ + ⌊y⌋ ≤ ⌊x + y⌋ ≤ ⌊x⌋ + ⌊y⌋ + 1 for all x and y.
Wikipedia — Kummer's theorem. For n ≥ m ≥ 0 and a prime p, the exponent of p in C(n, m) is the number of carries when m is added to n − m in base p, and equals (Sₚ(m) + Sₚ(n − m) − Sₚ(n))/(p − 1) with base-p digit sums Sₚ. Its example adds 3 = 11₂ and 7 = 111₂ with three carries for C(10, 3) = 120.
Wikipedia — Sieve of Eratosthenes. Crossing out may start at the square of each prime and stop once p² exceeds n; the sieve takes O(n log log n) operations, because the sum of the reciprocals of the primes grows like log log n, and the basic algorithm needs O(n) memory.
C++ — Euler231.cpp. The linked immutable revision contains parse_arguments, primes_up_to, exponent_in_factorial, solve and run_checkpoints; it is the implementation record for this article. Python — numeric types: integers have unlimited precision.

Problem 231 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int n = 20000000;
    int k = 15000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }

    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n) ||
            parse_int_after_prefix(arg, "--k=", options.k)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1 && options.k >= 0 && options.k <= options.n;
}

std::vector<int> primes_up_to(const int n) {
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(n + 1), 1);
    is_prime[0] = 0;
    is_prime[1] = 0;

    for (int i = 2; i * i <= n; ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            is_prime[static_cast<std::size_t>(j)] = 0;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

u64 exponent_in_factorial(const int n, const int p) {
    u64 exp = 0;
    int x = n;
    while (x > 0) {
        x /= p;
        exp += static_cast<u64>(x);
    }
    return exp;
}

u64 solve(const int n, const int k) {
    const int r = n - k;
    const std::vector<int> primes = primes_up_to(n);

    u64 sum = 0;
    for (int p : primes) {
        const u64 exp = exponent_in_factorial(n, p) -
                        exponent_in_factorial(k, p) -
                        exponent_in_factorial(r, p);
        sum += static_cast<u64>(p) * exp;
    }

    return sum;
}

bool run_checkpoints() {
    if (solve(10, 3) != 14ULL) {
        std::cerr << "Checkpoint failed for C(10,3)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.n, options.k) << '\n';
    return 0;
}

Python

def primes_up_to(n):
    is_prime = bytearray([1]) * (n + 1)
    is_prime[0] = is_prime[1] = 0
    
    for i in range(2, int(n**0.5) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = 0
                
    return [i for i in range(2, n + 1) if is_prime[i]]

def exponent_in_factorial(n, p):
    exp = 0
    x = n
    while x > 0:
        x //= p
        exp += x
    return exp

def solve(n=20000000, k=15000000):
    r = n - k
    primes = primes_up_to(n)
    
    total_sum = 0
    for p in primes:
        exp = exponent_in_factorial(n, p) - exponent_in_factorial(k, p) - exponent_in_factorial(r, p)
        total_sum += p * exp
        
    return str(total_sum)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;

public class Euler231 {
    static int[] primesUpTo(int n) {
        byte[] isPrime = new byte[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = 1;

        for (int i = 2; i * i <= n; i++) {
            if (isPrime[i] == 1) {
                for (int j = i * i; j <= n; j += i) {
                    isPrime[j] = 0;
                }
            }
        }

        int count = 0;
        for (int i = 2; i <= n; i++) {
            if (isPrime[i] == 1)
                count++;
        }

        int[] primes = new int[count];
        int idx = 0;
        for (int i = 2; i <= n; i++) {
            if (isPrime[i] == 1)
                primes[idx++] = i;
        }

        return primes;
    }

    static long exponentInFactorial(int n, int p) {
        long exp = 0;
        int x = n;
        while (x > 0) {
            x /= p;
            exp += x;
        }
        return exp;
    }

    public static String solve() {
        int n = 20000000;
        int k = 15000000;
        int r = n - k;

        int[] primes = primesUpTo(n);
        long sum = 0;

        for (int p : primes) {
            long exp = exponentInFactorial(n, p) - exponentInFactorial(k, p) - exponentInFactorial(r, p);
            sum += (long) p * exp;
        }

        return String.valueOf(sum);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}