Lucas' correspondence theorem, here or here, states that a binomial coefficient `Binomial[n,k]`

is equivalent to, mod prime `p`

, the product of binomial coefficients `Binomial[ni,ki]`

, where `ni`

and `ki`

are the digits in the base `p`

expansion of `n`

and `k`

, respectively.

The function `BinomialMod[n,k,p]`

does this calculation for **prime** modulus `p`

.

```
BinomialMod[n_, k_, p_] :=
With[{d = Max[IntegerLength[n, p], IntegerLength[k, p]]},
Mod[Apply[Times,
Mod[
Thread[Binomial[IntegerDigits[n, p, d], IntegerDigits[k, p, d]]],
p]],
p]]
```

As an example timing,

```
{AbsoluteTiming[Mod[Binomial[10^10, 10^5], 1123]],
AbsoluteTiming[BinomialMod[10^10, 10^5, 1123]]}
(* {{0.149513, 931}, {0.000054, 931}} *)
```

For a modulus `m`

equal to a product of **distinct** primes, the Chinese Remainder theorem gives the correct answer. The following is a brute force result.

```
AbsoluteTiming[Mod[Binomial[10^10, 10^5], 1009*1123*1259]]
(* {0.151880, 1131553699} *)
```

The fast way is with `ChineseRemainder`

.

```
With[{p = {1009, 1123, 1259}}, AbsoluteTiming[
ChineseRemainder[Map[BinomialMod[10^10, 10^5, #]&, p], p]]]]
(* {0.000281, 1131553699} *)
```

I hope someone can extend this answer to cases of a general composite modulus $m=p_1^{e_1} p_2^{e_2}...$ with repeated prime factors.

Great to see an answer to this after all this time :) It seems like Mathematica has nothing like this built in, so one would have to implement Granville's theorem or a different algorithm to the same end manually – Niklas B. – 2015-08-20T10:02:45.870