I need to perform some computations in a modulo ring, like

Mod[Subfactorial[n], m]
Mod[Binomial[n, k], m]

However, this is obviously much too slow for large n and k. Does Mathematica offer a possibility to perform those computations faster by using more efficient algorithms to compute the modulo? For example, there's Andrew Granville's generalization of Lucas' theorem, which would compute the binomial coefficient much faster.

I don't like to implement these by myself, because I want to use Mathematica to cross-check the results of another program of mine.


1 Answer 1


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]]}, 
             Thread[Binomial[IntegerDigits[n, p, d], IntegerDigits[k, p, d]]],

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.

  • $\begingroup$ 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 $\endgroup$
    – Niklas B.
    Aug 20, 2015 at 10:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.