Binomial coefficients modulo m

Today I (re)learned that you can compute binomial coefficients C(n,k) = n!/(k! (n-k)!), modulo some integer m, very efficiently. 

NOT as chatGPT or brave.AI or some stackoverflow people suggest, by computing n! and the denominator (mod m), and then multiply the former by the modular inverse of the latter:

That won't work because you will most likely (and certainly for n>m) get 0 for numerator and probably also the denominator.

However, when m is prime, we have the wonderful Lucas's theorem that tells us that 

C(n,k) = ∏ (C(n[i], k[i] ) ; i=1..L)  (mod m)       (where ∏ = product)

where n[i] and k[i] are the i-th digits of n and k in base p, and L is the length, i.e., number of base-p digits of k.

(For example, { np = digits(n,p), kp = digits(k,p) } in the PARI/GP language.)

Remark : In books and Wikipedia, you see the formula with "leading zeros" for k as to get the same length (i.e., L = max{ i | n[i] ≠ 0 }, and k[i] = 0 for the i beyond the "length" of k.
But given that C(n,0) = 1, those don't contribute and we can simply ignore the additional digits of n. 

Remark 2: Of course, we also assume n ≥ k, since C(n,k) = 0 for (integer) n < k. (For non-integer k that's an interesting, entirely different story, to be told another day...)

 

Now that we know C(n,k) (mod p) for primes p, we can also easily compute it for any squarefree modulus m. Indeed, this is exactly what does the Chinese Remainder Theorem: For given y[i] and m[i], this "theorem" (and algorithm) gives us x such that x = y[i] mod m[i] for all i.

Thus, in PARI/GP :

C_mod(n, k, m) = { if ( isprime(m),

   my ( kp = digits(k,m), np = digits(n,m)[-#kp..-1] );

   prod ( i=1,#kp, binomial(np[i],kp[i]), Mod(1,m) ), /* else: */

   my(f=factor(m));

   if ( vecmax(f[,2]) > 1, /* not squarefree */ C(n,k) % m,

        chinese([C_mod(n,k,p) | p <- f[,1]])

   )

)}

-- That's all! :-)

And this will be blazingly fast and use virtually no memory [for squarefree m], even for huge n & k where you won't even have enough RAM (and time) to compute the normal binomial(n,k).

PS: There is a way for the non-squarefree case, but that becomes more complicated -- see MaxAl's GP scripts here.