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/

fermat_binom.py

#!/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

fermat_binom_advanced.py

#!/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

lucas_binom.py

#!/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

lucas_crt_binom.py

#!/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

pascal_triangle.py

#!/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

pascal_triangle_row.py

#!/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

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).

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 👍