Computing Large Binomial Coefficients Modulo Prime / Non-Prime

Published on Jun 25, 2015

Last modified on Jun 28, 2015

Not rarely, in combinatoric problems it comes down to calculating the binomial coefficient

(\binom{n}{k})

for very large

n

and/or

k

modulo a number

m

. In general, the binomial coefficient can be formulated with factorials as

(\binom{n}{k}) = \frac{n!}{k! (n - k)!}, 0 \leq k \leq n

. The problem here is that factorials grow extremely fast which makes this formula computationally unsuitable because of quick overflows. For that reason, many problems in that category require the calculation of

(\binom{n}{k}) \mod m

. In this post I want to discuss ways to calculate the binomial coefficients for cases in which

m

is prime and when

m

is non-prime.

1. First simple approaches for any $m$

The good news is that there are easy ways to compute the binomial coefficient for any modulo

m

- the bad news is that they are not feasible for very large numbers. However, we will see later that we can break down very large numbers into sub-problems of much smaller numbers, in which case we can re-use techniques that we discuss now.

Using Pascal’s Triangle

The easiest (but also most limited) way to tackle this problem is to use a bottom-up dynamic programming approach by pre-computing Pascal’s Triangle. Pascal’s Triangle can be seen as a lookup table, where

(\binom{n}{k})

is the value at the

n

th row in the

k

th column. Pre-computing Pascal’s Triangle is very easy, since the current value is simply the sum of two neighboring cells from the previous row. We know that

(a + b) \mod m =

(a \mod m) +

(b \mod m) \mod m

, thus we can calculate each cell of Pascal’s Triangle with the modulo to avoid overflows.
Code: Computing Pascal’s Triangle, for any $m$

#!/usr/bin/env python3

"""

Calculating Pascal's Triangle to determine Binomial Coefficients nCk % m

"""

# Calculating Pascal's Triangle until row 'n' with modulo 'mod'

# bottom-up dynamic programming approach

# runtime: O(n**2)

# space: O(n**2)

def pascal_triangle(n,mod):

# prepare first row in table

pascal = []

pascal.append([1,0])

for i in range(1,n+1):

# initialize new row in table

pascal.append([])

pascal[i] = [0]*(i+2)

pascal[i][0] = 1

for j in range(1,i+1):

            # calculate current value based on neighbor values from previous row in triangle

pascal[i][j] = (pascal[i-1][j] + pascal[i-1][j-1])%mod

return pascal

if __name__ == '__main__':

n = 1009 # number of rows in pascal table

mod = 123456 # not a prime

# calculate pascal table with mod 123456

pascal = pascal_triangle(n,mod)

# (950 choose 100) mod 123456

print(pascal[950][100]) # should be 24942

view raw pascal_triangle.py hosted with ❤ by GitHub

Computing the triangle takes quadratic time

O (n^{2})

, since we compute each possible

k

for each row (

n

rows in total). Once the triangle is computed we can simply lookup the binomial coefficient in constant time

O (1)

for different

n

and

k

. The big problem arises when looking at the space complexity of

O (n^{2})

. For some problems

n

might be small enough, so there won’t be an issue, but for very large

n

we will not have enough memory to store the table, which means we have to find a better approach for those cases.
There is an easy upgrade in terms of space complexity. Each row from Pascal’s Triangle can be calculated by only knowing the previous row. This means we only have to keep track of the current row, thus the space complexity is linear

O (n)

.
Code: Computing single row in Pascal’s Triangle, for any $m$

#!/usr/bin/env python3

"""

Calculating a row in Pascal's Triangle to determine Binomial Coefficients

"""

# Calculating row 'n' of Pascal's Triangle with modulo 'mod' only remembering last calculated row

# bottom-up dynamic programming approach

# runtime: O(n**2)

# space: O(n)

def pascal_triangle_row(n,mod):

# initialize necessary space for row 'n'

pascal = [0]*(n+1)

pascal[0] = 1

for i in range(1,n+1):

# initialize auxiliary array

tmp = [0]*(n+1)

tmp[0] = 1

for j in range(1,i+1):

            # calculate current value based on neighbor values from previous row in triangle

tmp[j] = (pascal[j] + pascal[j-1])%mod

pascal = tmp

return pascal

if __name__ == '__main__':

n = 950 # row of pascal table

mod = 123456 # not a prime

# calculate pascal table row 950 with mod 123456

pascal_row = pascal_triangle_row(n,mod)

# (950 choose 100) mod 123456

print(pascal_row[100]) # should be 24942

view raw pascal_triangle_row.py hosted with ❤ by GitHub

Even though this approach allows to calculate larger binomial coefficients due to better space complexity, it still comes with a downside. In the first approach we can calculate the binomial coefficient in

O (1)

, after we have pre-computed Pascal’s Triangle in

O (n^{2})

. However in this approach, we need to recalculate the triangle or compute further rows each time a smaller or bigger row

n

is asked than the currently computed one. Consequently, this approach is expensive when we expect to calculate many different binomial coefficients that are not in the same row of Pascal’s Triangle.

2. Advanced approaches when $m$

is prime

In most related problems we are asked to calculate

(\binom{n}{k}) \mod m

, with

m

being a prime number. The choice of

m

being a prime number is by far no coincidence, since prime number modulus create a field in which certain Theorems hold.

Using Fermat’s Little Theorem

We remember that the binomial coefficient can be written as

\frac{n!}{k! (n - k)!} = \frac{n \cdot (n - 1) \cdot \dots \cdot (n - k + 1)}{k!}

. We can easily calculate the numerator of that equation, since

((a \mod m) \cdot (b \mod m))

\mod m = (a \cdot b) \mod m

. In order to divide the numerator by the denominator, we have to multiply the numerator with the multiplicative inverse of the denominator. Lucky for us, if

m

is prime, we can easily apply Fermat’s little theorem. Through this theorem we know that

a^{- 1} \equiv a^{p - 2} \mod p

, with

p

being a prime number. Consequently, we can easily find the multiplicative inverse through modular exponentiation of

(k!)^{m - 2}

.
Code: Using Fermat’s Theorem, for $m$
prime and $k < m$

#!/usr/bin/env python3

"""

Using Fermat's little theorem to calculate nCk mod m, for k < m and m is prime

Two versions:

1. Pre-Compute factorials and multiplicative inverses in O(n*logn) --> later lookup in O(1)

2. Compute directly --> no lookup --> each time O(n)

"""

# modular exponentiation: b^e % mod

# python's built-in pow(b,e,mod) would be faster than this function,

# I am still keeping it here for completion

def mod_exp(b,e,mod):

r = 1

while e > 0:

if (e&1) == 1:

r = (r*b)%mod

b = (b*b)%mod

e >>= 1

return r

# Using Fermat's little theorem to compute nCk mod p

# Note: p must be prime and k < p

def fermat_binom(n,k,p):

if k > n:

return 0

# calculate numerator

num = 1

for i in range(n,n-k,-1):

num = (num*i)%p

# calculate denominator

denom = 1

for i in range(1,k+1):

denom = (denom*i)%p

# numerator * denominator^(p-2) (mod p)

return (num * mod_exp(denom,p-2,p))%p

# Using Fermat's little theorem to pre-compute factorials and inverses

# Note: only works when p is prime and k < p

def fermat_compute(n,p):

facts = [0]*n

invfacts = [0]*n

facts[0] = 1

invfacts[0] = 1

for i in range(1,n):

# calculate factorial and corresponding inverse

facts[i] = (facts[i-1]*i)%p

invfacts[i] = mod_exp(facts[i],p-2,p)

return facts, invfacts

# Compute binomial coefficient from given pre-computed factorials and inverses

def binom_pre_computed(facts, invfacts, n, k, p):

# n! / (k!^(p-2) * (n-k)!^(p-2)) (mod p)

return (facts[n] * ((invfacts[k]*invfacts[n-k] % p))) % p

if __name__ == '__main__':

n = 1009 # number factorials to pre-compute

mod = 1000000007 # prime

# pre-compute factorials and inverses

facts, invfacts = fermat_compute(n,mod)

# print (950 choose 100) mod 1000000007 (with pre-computing)

print(binom_pre_computed(facts,invfacts,950,100,mod)) # should be 640644226

# (950 choose 100) mod 1000000007 (without pre-computing)

print(fermat_binom(950,100,mod)) # should be 640644226

view raw fermat_binom.py hosted with ❤ by GitHub

Computing the factorial

n!

takes

O (n)

and computing one multiplicative inverse takes

O (\log_{2} n)

. Note, that we implemented two versions here:
In the first version, we pre-computed the factorials and the corresponding inverses in an array at the cost of

O (n)

memory. This pre-computation takes

O (n \log_{2} n)

time and gives us the advantage of further lookups of binomial coefficients within the same range at constant time

O (1)

.
In the second version, without pre-computed lookups, we re-compute the values each time and only calculate the multiplicative inverse of one denominator, thus having a linear time

O (n)

in total each time we want to calculate a binomial coefficient.
Whether we use pre-computed tables or not highly depends on

n

’s expected size. In case we have very large

n

, pre-computing the factorial and inverse tables with Fermat’s theorem (space complexity

O (n)

and runtime of pre-computation

O (n \log_{2} n)

) might still be too expensive. In all other cases, pre-computing is desired since it will save us a lot of time later on. In the following I will only show algorithms without pre-computation, since those are also working for very large

n

, but changing those algorithms to pre-computed versions can be done easily (analogous to the current example), since the core idea is still the same.

Note, that this code is only working as long as

k < m

, since otherwise

k! \mod m = 0

, which means there is no inverse defined. Luckily, with slight modifications we can work around this issue. First we compare the degree of

m

in numerator and denominator, because the degree in the numerator has to be smaller or equal. Further, we cancel out any occurrence of

m

in both, numerator and denominator.
Code: Using Fermat’s Theorem, for $m$
prime

#!/usr/bin/env python3

"""

Using Fermat's little theorem to calculate nCk mod m, m is prime

Computation O(n)

"""

# modular exponentiation: b^e % mod

# python's built-in pow(b,e,mod) would be faster than this function,

# but I am still keeping it here for completion

def mod_exp(b,e,mod):

r = 1

while e > 0:

if (e&1) == 1:

r = (r*b)%mod

b = (b*b)%mod

e >>= 1

return r

# get degree of p in n! (exponent of p in the factorization of n!)

def fact_exp(n,p):

e = 0

u = p

t = n

while u <= t:

e += t//u

u *= p

return e

# Using Fermat's little theorem to compute nCk mod p

# considering cancelation of p in numerator and denominator

# Note: p must be prime

def fermat_binom_advanced(n,k,p):

# check if degrees work out

num_degree = fact_exp(n,p) - fact_exp(n-k,p)

den_degree = fact_exp(k,p)

if num_degree > den_degree:

return 0

if k > n:

return 0

# calculate numerator and cancel out occurrences of p

num = 1

for i in range(n,n-k,-1):

cur = i

while cur%p == 0:

cur //= p

num = (num*cur)%p

# calculate denominator and cancel out occurrences of p

denom = 1

for i in range(1,k+1):

cur = i

while cur%p == 0:

cur //= p

denom = (denom*cur)%p

# numerator * denominator^(p-2) (mod p)

return (num * mod_exp(denom,p-2,p))%p

if __name__ == '__main__':

mod = 1000000007 # prime

# (950 choose 100) mod 1000000007

print(fermat_binom_advanced(950,100,mod)) # should be 640644226

view raw fermat_binom_advanced.py hosted with ❤ by GitHub

We are now able to determine the binomial coefficient in linear time

O (n)

as long as the modulo is a prime number. We will see in the next sections how we can enhance the computation in terms of time complexity in some scenarios. Further, we will see how we can determine the binomial coefficients in case the modulo is not a prime number. The algorithms involving Fermat and Pascal are very crucial, since all further approaches that we are going to discuss in this post are relying on them to solve the core computation.

Using Lucas’ Theorem

Lucas’ theorem gives us a great way to break down the computation of one large binomial coefficient into the computation of smaller binomial coefficients. In order to do so, we first have to convert the digits of

n

and

k

into their base

m

representation

n = n_{0} + n_{1} \cdot m^{1} + n_{2} \cdot m^{2} \dots n_{x} \cdot m^{x}

and

k = k_{0} + k_{1} \cdot m^{1} + k_{2} \cdot m^{2} \dots k_{x} \cdot m^{x}

. We can then define sub problems:

(\binom{n}{k}) = \prod_{i = 0}^{x} (\binom{n_{i}}{k_{i}}) (\mod m)

. These problems can then be solved using Fermat’s theorem like we did before.
Code: Using Lucas’ + Fermat’s Theorem, for $m$
prime

#!/usr/bin/env python3

"""

Using Lucas' and Fermat's little theorem to calculate nCk mod m, m is prime

"""

# modular exponentiation: b^e % mod

# python's built-in pow(b,e,mod) would be faster than this function,

# but I am still keeping it here for completion

def mod_exp(b,e,mod):

r = 1

while e > 0:

if (e&1) == 1:

r = (r*b)%mod

b = (b*b)%mod

e >>= 1

return r

# get degree of p in n! (exponent of p in the factorization of n!)

def fact_exp(n,p):

e = 0

u = p

t = n

while u <= t:

e += t//u

u *= p

return e

# convert given number n into array of its base b representation

# most significant digit is at rightmost position in array

def get_base_digits(n,b):

d = []

while n > 0:

d.append(n % b)

n = n // b

return d

# Using Fermat's little theorem to compute nCk mod p

# considering cancelation of p in numerator and denominator

# Note: p must be prime

def fermat_binom_advanced(n,k,p):

# check if degrees work out

num_degree = fact_exp(n,p) - fact_exp(n-k,p)

den_degree = fact_exp(k,p)

if num_degree > den_degree:

return 0

if k > n:

return 0

# calculate numerator and cancel out occurrences of p

num = 1

for i in range(n,n-k,-1):

cur = i

while cur%p == 0:

cur //= p

num = (num*cur)%p

# calculate denominator and cancel out occurrences of p

denom = 1

for i in range(1,k+1):

cur = i

while cur%p == 0:

cur //= p

denom = (denom*cur)%p

# numerator * denominator^(p-2) (mod p)

return (num * mod_exp(denom,p-2,p))%p

# Using Lucas' theorem to split the problem into smaller sub-problems

# In this version assuming p is prime

def lucas_binom(n,k,p):

# get n and k in base p representation

np = get_base_digits(n,p)

kp = get_base_digits(k,p)

# calculate (nCk) = (n0 choose k0)*(n1 choose k1) ... (ni choose ki) (mod p)

binom = 1

for i in range(len(np)-1,-1,-1):

ni = np[i]

ki = 0

if i < len(kp):

ki = kp[i]

binom = (binom * fermat_binom_advanced(ni,ki,p)) % p

return binom

if __name__ == '__main__':

mod = 7 # prime

# (950 choose 100) mod 7

print(lucas_binom(950,100,mod)) # should be 2

view raw lucas_binom.py hosted with ❤ by GitHub

A little side fact: Lucas’ theorem gives an efficient way of determining whether our prime

m

is a divisor of the binomial coefficient, without calculating the coefficient itself. All we have to do is compare the digits of

n

and

k

in their base

m

representation

n = n_{0} + n_{1} \cdot m^{1} + n_{2} \cdot m^{2} \dots

and

k = k_{0} + k_{1} \cdot m^{1} + k_{2} \cdot m^{2} \dots

and if in at least one case

k_{i} > n_{i}

, then the product of all sub binomial coefficients would be 0, since

(\binom{n}{k}) = 0

, when

k > n

We used Lucas’ theorem to first split the problem into smaller binomial coefficients that we could then use to find our solution. In order to solve the small coefficients, we used Fermat’s theorem like in the previous section of this post. We could have also used Pascal’s Triangle in order to compute these small coefficients. Also, in case (n) and (k) are smaller than (m), Lucas’ theorem would not have brought any benefits, because the smallest sub problem would then be the original problem itself. However, in that case Lucas’ theorem does only bear the small overhead of bringing

n

and

k

in base

m

representation, which is why in general it is a good approach to first try splitting the problem into smaller sub problems with Lucas’ theorem, then using Pascal’s Triangle or Fermat’s theorem to solve these sub problems.

Remember, in this example we did not use pre-computation, which makes it more efficient when computing only one single very large

n

. Whether we use pre-computed tables or not highly depends on the expected size of

n

(which when using Lucas’ theorem, is the biggest prime factor in

m

) and whether we expect to calculate many different coefficients. In case our sub problems have still very large

n

, pre-computing the factorial and inverse tables with Fermat’s theorem (space complexity

O (n)

and initial runtime

O (n \log_{2} n)

) might still be too expensive to be feasible. In all other cases, pre-computing is desired since it will save us a lot of time later on.

All in all, now we have efficient ways of computing large binomial coefficients modulo a prime number.

3. $m$

’s prime factorization is square-free

What do we do in cases in which

m

turns out not to be a prime? The answer lies in the prime factorization of

m

. The following method works when

m

’s prime factorization is square-free, e.g.

6 = 2 \cdot 3

, which means

6

is square-free, whereas

24 = 2^{3} \cdot 3

is not square-free, since the prime

2

occurs more than once.
Lets assume

m

’s prime factors are

m = p_{0} \cdot p_{1} \dots \cdot p_{r}

. For each prime factor

p_{i}

m

, we can solve

(\binom{n}{k}) \mod p_{i}

, using the approaches we have discussed so far, which leaves us with

r

results

a_{0}, a_{1}, \dots, a_{r}

. Through the Chinese Remainder Theorem (CRT) we know that we can find a number

x

, which solves the following congruences

x \equiv a_{0} \mod p_{0}

x \equiv a_{1} \mod p_{1}

\dots

x \equiv a_{r} \mod p_{r}

.
Let’s define

m_{i} = \frac{p_{0} \cdot p_{1} \cdot \dots \cdot p_{r}}{p_{i}}

. We know that

m_{i}

and

p_{i}

must be coprime, then

g c d (m_{i}, p_{i}) = 1

, meaning

1 = s_{i} \cdot m_{i} + t_{i} \cdot p_{i}

. Using the Extended Euclidean Algorithm, we can find such

s_{i}

and

t_{i}

. Finally, we can find the solution that combines all the congruences to

x = \sum_{i = 0}^{r} (a_{i} \cdot m_{i} \cdot s_{i}) (\mod \prod_{i = 0}^{r} m_{i})

,
The CRT works as long as

p_{0} \dots p_{r}

are co-prime. We see, that by applying the CRT to those congruence systems, we can combine them into a single congruence modulo

m

, which solves our problem.
Code: Using Lucas’ + Fermat’s + Chinese Remainder Theorem, for (m) with square-free prime factors

#!/usr/bin/env python3

"""

Using Lucas' and Fermat's little theorem to calculate nCk mod m, m's prime factorization is square-free

Also using Chinese Remainder Theorem to combine congruences of prime factors.

"""

# modular exponentiation: b^e % mod

# python's built-in pow(b,e,mod) would be faster than this function,

# but I am still keeping it here for completion

def mod_exp(b,e,mod):

r = 1

while e > 0:

if (e&1) == 1:

r = (r*b)%mod

b = (b*b)%mod

e >>= 1

return r

# get degree of p in n! (exponent of p in the factorization of n!)

def fact_exp(n,p):

e = 0

u = p

t = n

while u <= t:

e += t//u

u *= p

return e

# convert given number n into array of its base b representation

# most significant digit is at rightmost position in array

def get_base_digits(n,b):

d = []

while n > 0:

d.append(n % b)

n = n // b

return d

# Extended Euclidean GCD

# compute x,y for ax + by = gcd(a,b)

# here, a,b are coprime, meaning gcd(a,b) = 1

def egcd(a, b):

x,y, u,v = 0,1, 1,0

while a != 0:

q, r = b//a, b%a

m, n = x-u*q, y-v*q

b,a, x,y, u,v = a,r, u,v, m,n

gcd = b

return (x, y)

# Chinese Remainder Theorem

# Combine given congruences to one solution

def crt(congruences):

# calculate the original modulo m

m = 1

for congruence in congruences:

m *= congruence[1]

# combine congruences

result = 0

for congruence in congruences:

s, t = egcd(m//congruence[1],congruence[1])

result += (congruence[0]*s*m)//congruence[1]

return result%m

# Using Fermat's little theorem to compute nCk mod p

# considering cancelation of p in numerator and denominator

# Note: p must be prime

def fermat_binom_advanced(n,k,p):

# check if degrees work out

num_degree = fact_exp(n,p) - fact_exp(n-k,p)

den_degree = fact_exp(k,p)

if num_degree > den_degree:

return 0

if k > n:

return 0

# calculate numerator and cancel out occurrences of p

num = 1

for i in range(n,n-k,-1):

cur = i

while cur%p == 0:

cur //= p

num = (num*cur)%p

# calculate denominator and cancel out occurrences of p

denom = 1

for i in range(1,k+1):

cur = i

while cur%p == 0:

cur //= p

denom = (denom*cur)%p

# numerator * denominator^(p-2) (mod p)

return (num * mod_exp(denom,p-2,p))%p

# Using Lucas' theorem to split the problem into smaller sub-problems

# p must be prime

def lucas_binom(n,k,p):

# get n and k in base p representation

np = get_base_digits(n,p)

kp = get_base_digits(k,p)

# calculate (nCk) = (n0 choose k0)*(n1 choose k1) ... (ni choose ki) (mod p)

binom = 1

for i in range(len(np)-1,-1,-1):

ni = np[i]

ki = 0

if i < len(kp):

ki = kp[i]

binom = (binom * fermat_binom_advanced(ni,ki,p)) % p

return binom

# Compute n choose k for given prime factors of m

# prime factors need to have multiplicity 1 in m

def binom(n,k,mod_facts):

# build congruences for all prime factors

congruences = []

for p in mod_facts:

# add (binom,p) to congruence list

congruences.append((lucas_binom(n,k,p),p))

# use CRT to combine congruences to one solution

return crt(congruences)

if __name__ == '__main__':

mod_facts = [3,5,7,11] # prime factors of m = 1155

# (8100 choose 4000) mod 1155

print(binom(8100,4000,mod_facts)) # should be 924

view raw lucas_crt_binom.py hosted with ❤ by GitHub

Finding the prime factorization of

m

is non-trivial and no efficient algorithm is known so far that can solve this problem for large

m

. However, as long as

m

has a reasonable size (and that is very likely in those combinatoric problems), finding the prime factorization is feasible. You can write your own algorithm or use some online calculators to determine the prime factors. Also remember that our current approach is only valid for

m

whose prime factors are square-free.

Now we are also able to compute large binomial coefficients modulo numbers whose prime factorization is square-free.

4. $m$

has prime factors with multiplicity higher than 1

Finally, lets have a look at all the remaining possible modulo numbers, namely those containing prime factors with a multiplicity higher than 1. For that purpose we can use Andrew Granville’s generalization of Lucas’ Theorem. This theorem also allows prime powers in the factorization of

m

. After applying the theorem, we could then use our other methods to calculate the sub problems and then like before use the Chinese Remainder Theorem to combine the sub solutions into one solution. With prime powers allowed, we can then factorize any modulus

m

, meaning we can find the binomial coefficients modulo any number

m

. I don’t have an implementation ready yet for the generalization of Lucas’ Theorem, but once I have some time I will try to add it to this post.

5. Quick Summary

We should do the following steps in order to compute large binomial coefficients

(\binom{n}{k}) \mod m

Find prime factors (and multiplicities) $p_{0}^{e_{0}}, \dots p_{r}^{e_{r}}$

m

with Fermat’s little theorem or Pascal’s Triangle
Use Chinese Remainder Theorem to combine sub results

Remember, if

p_{i} > n

, then there are no sub problems and we basically just solve the problem at hand using Fermat’s little theorem or Pascal’s Triangle.

That’s it! Now we can compute the binomial coefficients for very large numbers for any modulo

m

I use disqus as a comment system. If you click on the following button, then the disqus comment system will load in your browser and you agree to the disqus privacy policy. To delete your data from disqus you can contact their support team directly.

Load disqus comments

Use (generalized) Lucas’ Theorem to find all sub problems for each

(\binom{n}{k}) \mod p_{i}^{e_{i}}

Solve sub problems

(\binom{n}{k}) \mod p_{i}^{e_{i}}

No pain, No gain !!!

Thứ Năm, 11 tháng 4, 2019