Modulo of factorial divided by factorial

Question

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.

There are conditions for this to be an integer. Obviously, a>= max(b1,b2,...bm). Do you have more details - for example if b1+b2...bm=a then you are sure it is an integer because it is the number of ways to mix m groups of size bi without permutation inside the groups. — Jean Valj, Jan 15 '22 at 23:19
I know that a!/(b1! b2! ... bm!) is an integer, but since the numbers might be very big, we cannot always compute the factorials — Farah Mind, Jan 15 '22 at 23:26
How big is p? If not too big, a solution would be to consider number of occurrences of each prime (less than p) in the numerator and denominator. No idea if that would be optimal, but overflow could be avoided easily. — myrtlecat, Jan 15 '22 at 23:56
@user3386109 since p **is** prime, you can perform division by calculating modular inverses. I'll expand this into an answer. — myrtlecat, Jan 16 '22 at 00:09
[Removing common factors leaves a multiplication](https://stackoverflow.com/a/28058811/3386109). Then the multiplication can be performed by noting that `(a * b) mod p = ((a mod p) * (b mod p)) mod p`. — user3386109, Jan 16 '22 at 00:10
@myrtlecat I'm interested in seeing whether using modular inverses are faster than using GCDs to remove the common factors. — user3386109, Jan 16 '22 at 00:14
To remove common factors you could also get the prime factorizations of ```a!```, ```b1!```, ```b2!``` etc. and then subtract the exponents of the primes of ```a!``` by the exponents of the primes of each ```b_i!```. — Oli L, Jan 16 '22 at 00:16
@FarahMind Saying that the numbers "can be very big" is not helpful. What are the actual limits on `a`, `b` and `p`? — user3386109, Jan 16 '22 at 00:17
@user3386109 a and b can be from 40 to 60 and p is greater than 10000 — Farah Mind, Jan 16 '22 at 00:21
@KellyBundy OP only claimed that `a!` and `b!` were bigger than `long long int`, which is true — myrtlecat, Jan 16 '22 at 00:43
@KellyBundy 60! is bigger than unsigned long long int, is there another type that can supporte thsi order of number? — Farah Mind, Jan 16 '22 at 01:35
Then you can precompute at compile time the prime factors for all possible values of a and b and do it like Oli L suggested. As you confirmed, the fraction is an integer. Then you can combine as many factors as still fit into `long long int` and do `modulo` with `p` before multiplying the next prime factors of the integer. — Sebastian, Jan 16 '22 at 01:36
I would go with the previous suggestion with storing prime factors. You have about 19 prime factors below 60. You would have a 2D array storing those for 1..60 and you can initialize it with the correct values (by calculating externally and pasting into the code or by using `constexpr` functions, which do the calculations normally, but are executed at compile time). A second 1D array stores the prime numbers. As soon as you know the prime factors of the fraction result, it is a simple loop multiplying and from time to time do modulo (%). — Sebastian, Jan 16 '22 at 01:58
I would start with this code https://stackoverflow.com/questions/19019252/create-n-element-constexpr-array-in-c11/34465458#34465458 then change the array to two dimensions and instead of `arr[i] = i` insert a prime factor calculation `arr[i][p] = c`, where `c` is how often `i` can be divided by `p`. That is calculated at compile time and can be reasonably inefficient. — Sebastian, Jan 16 '22 at 02:06
@myrtlecat Nah, that's not really how they phrased it. They could've, but additionally said very big. Java for example has a class that easily supports integers with millions of digits, and even that isn't called `VeryBigInteger` :-). To me, the word "very" implies something much bigger than 82 digits (at least millions) and performance being an issue. — Kelly Bundy, Jan 16 '22 at 03:02

myrtlecat · Answer 1 · 2022-01-16T03:15:18.400

5

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.

edited Jan 16 '22 at 03:15

answered Jan 16 '22 at 00:42

myrtlecat

2,156
12
16

This answer makes it sound like canceling factors is optional, but in general you must first cancel at least factors of p. The reason is hinted at in other parts of the answer. Suppose (m/n) is an integer; then so is (p*m)/(p*n). But if you don't cancel p then the denominator is not invertible and the numerator is zero. – President James K. Polk Jan 16 '22 at 20:15
Yes, that's a good point. – myrtlecat Jan 16 '22 at 20:22
Since OP says that modulus is a prime, `inv(x)` can be `pow(x, p-2, p)` by Fermat's little theorem (using Python's `pow` function). – hilberts_drinking_problem Jan 17 '22 at 07:53
the third parameter being the modulus – Sebastian Jan 17 '22 at 07:58

score 4 · Answer 2 · edited Jan 16 '22 at 14:04

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)

Modulo of factorial divided by factorial

2 Answers2

Multiplication and modular arithmetic

Division and modular arithmetic