Doing computations in a modulo ring

10

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.

Niklas B.

Posted 2012-10-29T00:20:02.330

Reputation: 201

Answers

6

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.

KennyColnago

Posted 2012-10-29T00:20:02.330

Reputation: 14 269

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