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

Computing Large Binomial Coefficients Modulo Prime / Non-Prime

Computing Large Binomial Coefficients Modulo Prime / Non-Prime

Not rarely, in combinatoric problems it comes down to calculating the binomial coefficient
for very large and/or modulo a number . In general, the binomial coefficient can be formulated with factorials as . 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 . In this post I want to discuss ways to calculate the binomial coefficients for cases in which is prime and when is non-prime.

1. First simple approaches for any

The good news is that there are easy ways to compute the binomial coefficient for any modulo
- 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
is the value at the th row in the 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 , thus we can calculate each cell of Pascal’s Triangle with the modulo to avoid overflows.
Code: Computing Pascal’s Triangle, for any

#!/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
, since we compute each possible for each row ( rows in total). Once the triangle is computed we can simply lookup the binomial coefficient in constant time for different and . The big problem arises when looking at the space complexity of . For some problems might be small enough, so there won’t be an issue, but for very large 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
.
Code: Computing single row in Pascal’s Triangle, for any

#!/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
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
, after we have pre-computed Pascal’s Triangle in . However in this approach, we need to recalculate the triangle or compute further rows each time a smaller or bigger row 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

is prime

In most related problems we are asked to calculate
, with being a prime number. The choice of 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
. We can easily calculate the numerator of that equation, since . 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 is prime, we can easily apply Fermat’s little theorem. Through this theorem we know that , with being a prime number. Consequently, we can easily find the multiplicative inverse through modular exponentiation of .
Code: Using Fermat’s Theorem, for
prime and
#!/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
takes and computing one multiplicative inverse takes . 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
memory. This pre-computation takes time and gives us the advantage of further lookups of binomial coefficients within the same range at constant time .
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
in total each time we want to calculate a binomial coefficient.
Whether we use pre-computed tables or not highly depends on
’s expected size. In case we have very large , pre-computing the factorial and inverse tables with Fermat’s theorem (space complexity and runtime of pre-computation ) 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
, 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
, since otherwise , which means there is no inverse defined. Luckily, with slight modifications we can work around this issue. First we compare the degree of in numerator and denominator, because the degree in the numerator has to be smaller or equal. Further, we cancel out any occurrence of in both, numerator and denominator.
Code: Using Fermat’s Theorem, for
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
We are now able to determine the binomial coefficient in linear time
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
and into their base representation and . We can then define sub problems: . These problems can then be solved using Fermat’s theorem like we did before.
Code: Using Lucas’ + Fermat’s Theorem, for
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
is a divisor of the binomial coefficient, without calculating the coefficient itself. All we have to do is compare the digits of and in their base representation and and if in at least one case , then the product of all sub binomial coefficients would be 0, since , when
.
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
and in base 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
. Whether we use pre-computed tables or not highly depends on the expected size of (which when using Lucas’ theorem, is the biggest prime factor in ) and whether we expect to calculate many different coefficients. In case our sub problems have still very large , pre-computing the factorial and inverse tables with Fermat’s theorem (space complexity and initial runtime
) 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.

’s prime factorization is square-free

What do we do in cases in which
turns out not to be a prime? The answer lies in the prime factorization of . The following method works when ’s prime factorization is square-free, e.g. , which means is square-free, whereas is not square-free, since the prime occurs more than once.
Lets assume
’s prime factors are . For each prime factor of , we can solve , using the approaches we have discussed so far, which leaves us with results . Through the Chinese Remainder Theorem (CRT) we know that we can find a number , which solves the following congruences , , , .
Let’s define
. We know that and must be coprime, then , meaning . Using the Extended Euclidean Algorithm, we can find such and . Finally, we can find the solution that combines all the congruences to ,
The CRT works as long as
are co-prime. We see, that by applying the CRT to those congruence systems, we can combine them into a single congruence modulo , 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
is non-trivial and no efficient algorithm is known so far that can solve this problem for large . However, as long as 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
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.

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
. 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 , meaning we can find the binomial coefficients modulo any number . 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
:
  • Find prime factors (and multiplicities)
of
  • with Fermat’s little theorem or Pascal’s Triangle
  • Use Chinese Remainder Theorem to combine sub results
Remember, if
, 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
.

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
  • Solve sub problems
  • Không có nhận xét nào:

    Đăng nhận xét

    Bài G - Educatioal Round 62

    Đề bài: Bạn được cho 1 đồ thị vô hướng đặc biệt. Nó bao gồm $2n$ đỉnh được đánh số từ 1 đến 2n. Dưới đây là một số đặc tính của đồ thị: + ...