Skip to content

Instantly share code, notes, and snippets.

@fishi0x01
Last active March 21, 2020 15:19
Show Gist options
  • Save fishi0x01/46ebb1d310df5f9b5d8c to your computer and use it in GitHub Desktop.
Save fishi0x01/46ebb1d310df5f9b5d8c to your computer and use it in GitHub Desktop.
Calculating large binomial coefficients modulo prime / non-prime numbers (nCk mod m). Code for blog post https://fishi.devtail.io/weblog/2015/06/25/computing-large-binomial-coefficients-modulo-prime-non-prime/
#!/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
#!/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
#!/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
#!/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
#!/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
#!/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
@ksenobojca
Copy link

Python has pow with modulo argument which does what mod_exp does. Plus it's compiled so will work faster (unless you use stuff like PyPy).

@fishi0x01
Copy link
Author

fishi0x01 commented Feb 10, 2019

Good point!
Will keep the mod_exp for completion, but with a comment annotated above which states what you mention.
Thank you for the hint 👍

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment