What is the most efficient way to compute factorials modulo a prime?

Question

Do you know any algorithm that calculates the factorial after modulus efficiently?

For example, I want to program:

for(i=0; i<5; i++)
  sum += factorial(p-i) % p;

But, p is a big number (prime) for applying factorial directly $(p \leq 10^ 8)$.

In Python, this task is really easy, but i really want to know how to optimize.

Answer 1

(This answer was initially posted by the asker jonaprieto inside the question.)

I remember Wilson's theorem, and I noticed little things:

In the above program, it is better if I write: $$\begin{align} (p-1)! &\equiv -1 &\pmod p\\ (p-2)! &\equiv (p-1)! (p-1)^ {-1} \equiv \bf{1} &\pmod p\\ (p-3)! &\equiv (p-2)! (p-2)^ {-1} \equiv \bf{(p-2)^{-1}} &\pmod p\\ (p-4)! &\equiv (p-3)! (p-3)^ {-1} \equiv \bf{(p-2)^{-1}} \bf{(p-3)^{-1}} &\pmod p\\ \ (p-5)! &\equiv (p-4)! (p-4)^ {-1} \equiv \bf{(p-2)^{-1}} \bf{(p-3)^{-1}}\bf{(p-4)^{-1}} &\pmod p\\ \end{align}$$

And you can find $(p-i)^{-1}$ because $\operatorname{gcd}(p, p-i) = 1$, so with the extended Euclidian algorithm you can find the value of $(p-i)^{-1}$, that is the inverse modulus.

You can view the same congruences too, like to: $$\begin{align*} (p-5)! &\equiv (p-24)^{-1}&\pmod p\\ (p-4)! &\equiv (p+6)^{-1}&\pmod p\\ (p-3)! &\equiv (p-2)^{-1} &\pmod p\\ (p-2)! &\equiv 1&\pmod p\\ (p-1)! &\equiv -1&\pmod p\\ \end{align*} $$ so, the sum is equal: $$ (-24)^{-1}+(6)^{-1} +(-2)^{-1}$$ and if you factorize in the beginning the factorials you get $$ 8\cdot (-24)^{-1} \pmod p$$ And, voila, inverse modulus is more efficient than factorials.

Answer 2

The example that you are posting is very closely related to Euler problem #381. So I will post an answer that doesn't solve the Euler problem. I will post how you can calculate factorials modulo a prime.

So: How to calculate n! modulo p?

Quick observation: If n ≥ p, then n! has a factor p, so the result is 0. Very quick. And if we ignore the requirement that p should be a prime then let q be the smallest prime factor of p, and n! modulo p is 0 if n ≥ q. There's also not much reason to require that p is a prime to answer your question.

Now in your example (n - i)! for 1 ≤ i ≤ 5 came up. You don't have to calculate five factorials: You calculate (n - 5)!, multiply by (n - 4) go get (n - 4)!, multiply by (n - 3) to get (n - 3)! etc. This reduces the work by almost a factor 5. Don't solve the problem literally.

The question is how to calculate n! modulo m. The obvious way is to calculate n!, a number with roughly n log n decimal digits, and calculate the remainder modulo p. That's hard work. Question: How can we get this result quicker? By not doing the obvious thing.

We know that ((a * b * c) modulo p = (((a * b) modulo p) * c) modulo p.

To calculate n!, we would normally start with x = 1, then multiply x by 1, 2, 3, ... n. Using the modulo formula, we calculate n! modulo p without calculating n!, by starting with x = 1, and then for i = 1, 2, 3, .., n we replace x with (x * i) modulo p.

We always have x < p and i < n, so we only need enough precision to calculate x * p, not the much higher precision to calculate n!. So to calculate n! modulo p for p ≥ 2 we take the following steps:

Step 1: Find the smallest prime factor q of p. If n ≥ q then the result is 0.
Step 2: Let x = 1, then for 1 ≤ i ≤ n replace x with (x * i) modulo p, and x is the result. 

(Some answers mention Wilson's theorem, which only answers the question in the very special case of the example given, and is very useful to solve Euler problem #381, but in general isn't useful to solve the question that was asked).

Answer 3

This is my implementation use of the wilson's theorem:

The factMOD function is the one to call to compute (n!) % MOD when MOD-n is little against n.

Do someone know an other efficient approach when it's not the case (e.g: n=1e6 and MOD=1e9+7) ?

ll powmod(ll a, ll b){//a^b % MOD
  ll x=1,y=a;
  while(b){
    if(b&1){
      x*=y; if(x>=MOD)x%=MOD;
    }
    y*=y; if(y>=MOD)y%=MOD;
    b>>=1;
  }
  return x;
} 
ll InverseEuler(ll n){//modular inverse of n
  return powmod(n,MOD-2);
}
ll factMOD(ll n){ //n! % MOD efficient when MOD-n<n
   ll res=1,i;
   for(i=1; i<MOD-n; i++){
     res*=i;
     if(res>=MOD)res%=MOD;
   }
   res=InverseEuler(res);   
    if(!(n&1))
      res= -res +MOD;
  }
  return res%MOD;
}