Modulo of factorial divided by factorial-CodePudding

How to calculate a!/(b1! b2! ... bm!) modulo p, where p is a prime number? The factorial of a and b can be very big (long long int is not sufficient) so I need to pass to modulo.

CodePudding user response：

If a, bs and p are fairly small, prefer @KellyBundy's approach of cancelling factors, or counting prime factors.

Multiplication and modular arithmetic

Given integers m and n and some other integer k:

(m * n) modulo k = ((m modulo k) * (n mod k)) modulo k

This allows a large product to be calculated modulo p without worrying about overflow, since we can always keep the arguments in the range [0, k).

For example to compute the factorial a! modulo k, in python:

def fact(a, k):
    if a == 0:
        return 1
    else:
        return ((a % k) * fact(a - 1, k)) % k

Division and modular arithmetic

If p is a prime then for any integer n that is not divisible by p, we can find an integer which I'll call inv(n) such that:

(n * inv(n)) modulo p = 1

This number is called the modular inverse of n. There are various algorithms to find modular inverses, which I won't describe here (but see e.g. here).

Now, given integers n and m, and assuming that m / n is an integer, we can apply the rule:

(m / n) modulo p = (m * inv(n)) modulo p

So provided we can calculate modular inverses, we can convert division to multiplication, and then apply the previous rule.

CodePudding user response：

Another way, listing the factors 1 to a, then canceling with all divisors, then multiplying modulo p:

#include <iostream>
#include <vector>

int gcd(int a, int b) {
  return b ? gcd(b, a % b) : a;
}

int main() {
  int a = 60;
  std::vector<int> bs = {13, 7, 19};
  int p = 10007;

  std::vector<int> factors(a);
  for (int i=0; i<a; i  )
    factors[i] = i   1;
  for (int b : bs) {
    while (b > 1) {
      int d = b--;
      for (int& f : factors) {
        int g = gcd(f, d);
        f /= g;
        d /= g;
      }
    }
  }
  int result = 1;
  for (int f : factors)
    result = result * f % p;
  std::cout << result;
}

Prints 5744, same as this Python code:

from math import factorial, prod

a = 60
bs = [13, 7, 19]
p = 10007

num = factorial(a)
den = prod(map(factorial, bs))
print(num // den % p)