Not rarely, in combinatoric problems it comes down to calculating the
binomial coefficient (nk)
for very large
n and/or
k modulo a number
m. In general, the binomial coefficient can be formulated with
factorials as
(nk)=n!k!(n−k)!,0≤k≤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
(nk)modm. 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
(nk)
is the value at the
nth row in the
kth
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)modm=(amodm)+(bmodm)modm, thus we can calculate each cell of Pascal’s Triangle with the modulo to avoid overflows.
Code: Computing Pascal’s Triangle, for any m
| 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 |
| def pascal_triangle(n,mod): |
| # prepare first row in table |
| # initialize new row in table |
| # 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 |
| 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 |
Computing the triangle takes quadratic time
O(n2)
, 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(n2). 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
| 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 |
| def pascal_triangle_row(n,mod): |
| # initialize necessary space for row 'n' |
| # initialize auxiliary array |
| # calculate current value based on neighbor values from previous row in triangle |
| tmp[j] = (pascal[j] + pascal[j-1])%mod |
| 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 |
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(n2). 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
(nk)modm
, 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
n!k!(n−k)!=n⋅(n−1)⋅⋯⋅(n−k+1)k!
. We can easily calculate the numerator of that equation, since
((amodm)⋅(bmodm))modm=(a⋅b)modm.
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≡ap−2modp, 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
| Using Fermat's little theorem to calculate nCk mod m, for k < m and m is prime |
| 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 |
| # Using Fermat's little theorem to compute nCk mod p |
| # Note: p must be prime and k < p |
| for i in range(n,n-k,-1): |
| # 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 |
| # calculate factorial and corresponding inverse |
| facts[i] = (facts[i-1]*i)%p |
| invfacts[i] = mod_exp(facts[i],p-2,p) |
| # 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 |
| # 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 |
Computing the factorial n!
takes
O(n) and computing one multiplicative inverse takes
O(log2n). 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(nlog2n) 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(nlog2n))
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!modm=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
| Using 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 |
| # get degree of p in n! (exponent of p in the factorization of n!) |
| # Using Fermat's little theorem to compute nCk mod p |
| # considering cancelation of p in numerator and denominator |
| 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: |
| # calculate numerator and cancel out occurrences of p |
| for i in range(n,n-k,-1): |
| # calculate denominator and cancel out occurrences of p |
| # numerator * denominator^(p-2) (mod p) |
| return (num * mod_exp(denom,p-2,p))%p |
| if __name__ == '__main__': |
| # (950 choose 100) mod 1000000007 |
| print(fermat_binom_advanced(950,100,mod)) # should be 640644226 |
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=n0+n1⋅m1+n2⋅m2…nx⋅mx and
k=k0+k1⋅m1+k2⋅m2…kx⋅mx. We can then define sub problems:
(nk)=∏xi=0(niki)(modm). These problems can then be solved using Fermat’s theorem like we did before.
Code: Using Lucas’ + Fermat’s Theorem, for m
prime
| 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 |
| # get degree of p in n! (exponent of p in the factorization of n!) |
| # 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): |
| # Using Fermat's little theorem to compute nCk mod p |
| # considering cancelation of p in numerator and denominator |
| 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: |
| # calculate numerator and cancel out occurrences of p |
| for i in range(n,n-k,-1): |
| # calculate denominator and cancel out occurrences of 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 |
| # 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) |
| for i in range(len(np)-1,-1,-1): |
| binom = (binom * fermat_binom_advanced(ni,ki,p)) % p |
| if __name__ == '__main__': |
| print(lucas_binom(950,100,mod)) # should be 2 |
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=n0+n1⋅m1+n2⋅m2… and
k=k0+k1⋅m1+k2⋅m2… and if in at least one case
ki>ni, then the product of all sub binomial coefficients would be 0, since
(nk)=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(nlog2n)
)
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⋅3, which means
6 is square-free, whereas
24=23⋅3 is not square-free, since the prime
2 occurs more than once.
Lets assume
m
’s prime factors are
m=p0⋅p1⋯⋅pr. For each prime factor
pi of
m, we can solve
(nk)modpi, using the approaches we have discussed so far, which leaves us with
r results
a0,a1,…,ar. Through the
Chinese Remainder Theorem (CRT) we know that we can find a number
x, which solves the following congruences
x≡a0modp0,
x≡a1modp1,
…,
x≡armodpr.
Let’s define
mi=p0⋅p1⋅⋯⋅prpi
. We know that
mi and
pi must be coprime, then
gcd(mi,pi)=1, meaning
1=si⋅mi+ti⋅pi. Using the
Extended Euclidean Algorithm, we can find such
si and
ti. Finally, we can find the solution that combines all the congruences to
x=∑ri=0(ai⋅mi⋅si)(mod∏ri=0mi),
The CRT works as long as
p0…pr
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
| 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 |
| # get degree of p in n! (exponent of p in the factorization of n!) |
| # 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): |
| # compute x,y for ax + by = gcd(a,b) |
| # here, a,b are coprime, meaning gcd(a,b) = 1 |
| b,a, x,y, u,v = a,r, u,v, m,n |
| # Chinese Remainder Theorem |
| # Combine given congruences to one solution |
| # calculate the original modulo m |
| for congruence in congruences: |
| for congruence in congruences: |
| s, t = egcd(m//congruence[1],congruence[1]) |
| result += (congruence[0]*s*m)//congruence[1] |
| # Using Fermat's little theorem to compute nCk mod p |
| # considering cancelation of p in numerator and denominator |
| 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: |
| # calculate numerator and cancel out occurrences of p |
| for i in range(n,n-k,-1): |
| # calculate denominator and cancel out occurrences of 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 |
| # 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) |
| for i in range(len(np)-1,-1,-1): |
| binom = (binom * fermat_binom_advanced(ni,ki,p)) % p |
| # 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 |
| # add (binom,p) to congruence list |
| congruences.append((lucas_binom(n,k,p),p)) |
| # use CRT to combine congruences to one solution |
| 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 |
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
(nk)modm
:
- Find prime factors (and multiplicities) pe00,…perr
of
m
- with Fermat’s little theorem or Pascal’s Triangle
- Use Chinese Remainder Theorem to combine sub results
Remember, if pi>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.
Use (generalized) Lucas’ Theorem to find all sub problems for each (nk)modpeii
Solve sub problems (nk)modpeii
