algorithm - 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).