Fast way to calculate n! mod m where m is prime?

Question

I was curious if there was a good way to do this. My current code is something like:

def factorialMod(n, modulus):     ans=1     for i in range(1,n+1):         ans = ans * i % modulus         return ans % modulus 

But it seems quite slow!

I also can't calculate n! and then apply the prime modulus because sometimes n is so large that n! is just not feasible to calculate explicitly.

I also came across http://en.wikipedia.org/wiki/Stirling%27s_approximation and wonder if this can be used at all here in some way?

Or, how might I create a recursive, memoized function in C++?

Answer 1

n can be arbitrarily large

Well, n can't be arbitrarily large - if n >= m, then n! ≡ 0 (mod m) (because m is one of the factors, by the definition of factorial).

---

Assuming n << m and you need an exact value, your algorithm can't get any faster, to my knowledge. However, if n > m/2, you can use the following identity (Wilson's theorem - Thanks @Daniel Fischer!)

to cap the number of multiplications at about m-n

 (m-1)! ≡ -1 (mod m) 1 * 2 * 3 * ... * (n-1) * n * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m) n! * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m) n! ≡ -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m) 

This gives us a simple way to calculate n! (mod m) in m-n-1 multiplications, plus a modular inverse:

 def factorialMod(n, modulus):     ans=1     if n <= modulus//2:         #calculate the factorial normally (right argument of range() is exclusive)         for i in range(1,n+1):             ans = (ans * i) % modulus        else:         #Fancypants method for large n         for i in range(n+1,modulus):             ans = (ans * i) % modulus         ans = modinv(ans, modulus)         ans = -1*ans + modulus     return ans % modulus 

---

We can rephrase the above equation in another way, that may or may-not perform slightly faster. Using the following identity:

we can rephrase the equation as

 n! ≡ -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m) n! ≡ -[(n+1-m) * ... * (m-2-m) * (m-1-m)]-1 (mod m)        (reverse order of terms) n! ≡ -[(-1) * (-2) * ... * -(m-n-2) * -(m-n-1)]-1 (mod m) n! ≡ -[(1) * (2) * ... * (m-n-2) * (m-n-1) * (-1)(m-n-1)]-1 (mod m) n! ≡ [(m-n-1)!]-1 * (-1)(m-n) (mod m) 

This can be written in Python as follows:

 def factorialMod(n, modulus):     ans=1     if n <= modulus//2:         #calculate the factorial normally (right argument of range() is exclusive)         for i in range(1,n+1):             ans = (ans * i) % modulus        else:         #Fancypants method for large n         for i in range(1,modulus-n):             ans = (ans * i) % modulus         ans = modinv(ans, modulus)          #Since m is an odd-prime, (-1)^(m-n) = -1 if n is even, +1 if n is odd         if n % 2 == 0:             ans = -1*ans + modulus     return ans % modulus 

---

If you don't need an exact value, life gets a bit easier - you can use Stirling's approximation to calculate an approximate value in O(log n) time (using exponentiation by squaring).

---

Finally, I should mention that if this is time-critical and you're using Python, try switching to C++. From personal experience, you should expect about an order-of-magnitude increase in speed or more, simply because this is exactly the sort of CPU-bound tight-loop that natively-compiled code excels at (also, for whatever reason, GMP seems much more finely-tuned than Python's Bignum).

Answer 2

Expanding my comment to an answer:

Yes, there are more efficient ways to do this. But they are extremely messy.

So unless you really need that extra performance, I don't suggest to try to implement these.

---

The key is to note that the modulus (which is essentially a division) is going to be the bottleneck operation. Fortunately, there are some very fast algorithms that allow you to perform modulus over the same number many times.

• Division by Invariant Integers using Multiplication
• Montgomery Reduction

These methods are fast because they essentially eliminate the modulus.

---

Those methods alone should give you a moderate speedup. To be truly efficient, you may need to unroll the loop to allow for better IPC:

Something like this:

ans0 = 1 ans1 = 1 for i in range(1,(n+1) / 2):     ans0 = ans0 * (2*i + 0) % modulus         ans1 = ans1 * (2*i + 1) % modulus      return ans0 * ans1 % modulus 

but taking into account for an odd # of iterations and combining it with one of the methods I linked to above.

Some may argue that loop-unrolling should be left to the compiler. I will counter-argue that compilers are currently not smart enough to unroll this particular loop. Have a closer look and you will see why.

---

Note that although my answer is language-agnostic, it is meant primarily for C or C++.