Modular Multiplicative Inverse and Modular Combinatorics (nCr % prime) in C++

In mathematics, particularly in the area of arithmetic, a modular multiplicative inverse of an integer a is an integer x such that the product ax is congruent to 1 with respect to the modulus m. Written as:

which is the shorthand way of writing the statement that the remainder after dividing ax by the integer m is 1. It is important to note that the modular inverse exists if and only if a and m are coprime (ie gcd(a, m) = 1).

Finding the modular inverse in O(m)

The problem is to compute the modular inverse modulo m for every number upto n . We denote the modular inverse of i by inv[i] . Then for i > 1 , when modulus m is prime, the following equation is valid:

inv[i] = ((m - \left\lfloor \frac{m}{i} \right\rfloor) \cdot inv[m \bmod i]) \bmod m

Implementation in C++

    inv[1] = 1;

    for (int i=2; i<= n; i++) {
        inv[i] = ((m - (m/i))*inv[m%i])%m;
    }

Derivation

Using Division Algorithm,

m = q \cdot i + r

where,

• q = \left\lfloor \frac{m}{i} \right\rfloor is the quotient
• r = m \bmod i is the remainder

Rearranging the terms, we get

r = m - q \cdot i

multiplying both sides by inv[r] (assuming inv[r] exists)

r \cdot inv[r] = m \cdot inv[r] - q \cdot i \cdot inv[r]

since r \cdot inv[r] \equiv 1 \, (\bmod \, m) , we get,

1 \equiv m \cdot inv[r] - q \cdot i \cdot inv[r] \, (\bmod \, m)

m \cdot inv[r] \equiv 0 \, (\bmod \, m) , this is because, m \bmod m is 0

hence, we get,

i \cdot (-q \cdot inv[r]) \equiv 1 \, (\bmod \, m)

So, i^{-1} can be written as:

i^{-1} = (-q) \cdot inv[r] \, \bmod \, m

substituting q and r ,

i^{-1} = (-\left\lfloor \frac{m}{i} \right\rfloor) \cdot inv[m \bmod i] \, \bmod \, m

Using the properties of modulo, we can write the positive equivalent mod m as

(-x) \equiv (m - x) \, (\bmod \, m)

so our equation becomes,

i^{-1} = (m - \left\lfloor \frac{m}{i} \right\rfloor) \cdot inv[m \bmod i] \, \bmod \, m

which is

inv[i] = ((m - (m/i)) * inv[m \% i]) \% m

Calculating Inverse Factorials Modulo Prime

We know that,

i! = (i - 1)! \cdot i

Taking the modular inverse of both sides,

(i!)^{-1} = ((i - 1)!)^{-1} \cdot (i)^{-1} \bmod p

Inverse factorial can now be calculated using this recurrence relation,

Computing ⁿC_r % prime

We know that:

nCr =  n! / ((r!) x (n-r)!) 
or
nCr =  factorial(n) / (factorial(r) x factorial(n-r))

Now, using modular inverse, we can write

nCr % p = (factorial[n]*modInverse(factorial[r]) % p *
               modInverse(factorial[n-r]) % p) % p;

Here modInverse() means modular inverse under modulo p.

This is helpful in computing nCr % prime when both n and r are large. Since computing n! and r! will cause Integer overflow.

Code Implementation

//
//  main.cpp
//  Modular Combinatorics
//
//  Created by Himanshu on 18/02/22.
//

#include <iostream>
#define p 1000000007
using namespace std;
typedef long long ll;

void precompute (int n, ll fact[], ll inv[], ll invFact[]) {
    fact[0] = 1;
    
    for (int i=1; i<= n; i++) {
        fact[i] = (fact[i-1]*i)%p;
    }
    
    inv[1] = 1;
    for (int i=2; i<= n; i++) {
        inv[i] = ((p - (p/i))*inv[p%i])%p;
    }
    
    invFact[0] = invFact[1] = 1;
       
    for (int i=2; i<= n; i++) {
        invFact[i] = (invFact[i-1]*inv[i])%p;
    }
}


ll computeCombinatorics (int n, int r, ll fact[], ll inv[], ll invFact[]) {
    ll ans = (((fact[n]*invFact[r])%p)*(invFact[n-r]%p))%p;
    return ans;
}

int main (int argc, const char * argv[]) {
    int n = 100000;
    int r = 500;
    
    ll *fact = new ll[n]();
    ll *inv = new ll[n]();
    ll *invFact = new ll[n]();
    
    precompute(n, fact, inv, invFact);
    cout<<computeCombinatorics(n, r, fact, inv, invFact)<<endl;
    
    return 0;
}

Output

115855765

Time Complexity: O(n)

Here’s a working example: Ideone

References:
Modulo Inverse
Compute nCr % p (Using Fermat Little Theorem)