Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question

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]]}, 
      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.

share|improve this answer
    
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. Aug 20 at 10:02

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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