How to solve nCr%p using Fermat’s Little Theorem?

Share your love

Problem statement

Given three numbers n, r and p, compute the value of nCr%p. Here p is a prime number greater than n, and nCr is the Binomial Coefficient.

Example: 

Input:  n = 10, r = 2, p = 13

Output: 6

Explanation: 10C2 is 45 and 45 % 13 is 6.

In the previous blog, we have seen how we can solve nCr%p with the help of Dynamic Programming. With this approach, we can solve nCr %p for both the cases – when p is a prime number and when p isn’t a prime number. 

The drawback of that solution was its time complexity was O(n*r) which cannot be optimized for solutions when p is a non-prime number. But when p is a prime number, it can be further optimized to O(n) using Fermat’s Little Theorem.

Fermat’s Little Theorem

Fermat’s little theorem states that if p is a prime number, then for any integer a, the number ap – a is an integer multiple of p. In the notation of modular arithmetic, this is expressed as: 

ap = a (mod p)

If a is not divisible by p, Fermat’s little theorem is equivalent to the statement a p – 1 – 1 is an integer multiple of p, then

ap-1 = 1 (mod p)

If we multiply both sides by a-1, we get. 

ap-2 = a-1 (mod p)

Approach

The idea is to replace all the values of the inverse mod in the nCr formula using the above formula of Fermat’s Little Theorem.

Thus the formula of ncr becomes

nCr = (factorial(n)%p * factorial(r)p-2%p * factorial(n-r)p-2%p)%p

Algorithm

  1. Initialize an array of n.
  2. iterate through it from the first index till n and find the factorial of each number. Fix the value of the 0th index as 1
  3. Use these factorials to find the value of the ncr using the above Fermat’s formula.
  4. Calculate the value of the inverse mod using the power function
  5. Check whether the power is even or odd. If odd, return the answer accordingly.

We have used bit manipulation to solve the power function as it will give us the answer in logarithmic time complexity as any other method will take linear time to process and we might get a stack overflow error and long instead of int to avoid overflow.

JAVA
package com.Tekolio.Extras;
// solve nCr%p
public class solve_cobnimatrics_1 {
    public static void main(String[] args) {
        int n = 41;
        int r = 27;
        int p = 139;
        System.out.print("Value of "+n+"C"+r+" % "+p+" is "+solve(n, r, p));
    }
public static int solve(int A, int B, int C) {
        //Base Case
        if (B == 0)
            return 1;
        //Precomputing the factorial values
        long[] fac = new long[A + 1];
        fac[0] = 1;

        for (int i = 1; i <= A; i++)
            fac[i] = fac[i - 1] * i % C;
        
        long ncr = ((fac[A] * modInverse(fac[B], C)
                % C * modInverse(fac[A - B], C)
                % C)% C);
        return (int)ncr;
    }
    static long modInverse(long n, long p)
    {
        return power(n, p - 2, p);
    }
    static long power(long x, long y, long p)
    {

        // Initialize result
        long res = 1;

        // Update x if it is more than or
        // equal to p
        x = x % p;

        while (y > 0) {

            // If y is odd, multiply x
            // with result
            if (y % 2 == 1)
                res = (res * x) % p;

            // y must be even now
            y = y >> 1; // y = y/2
            x = (x * x) % p;
        }

        return res;
    }
}
C++
#include <bits/stdc++.h>
using namespace std;
 
/* Iterative Function to calculate (x^y)%p
in O(log y) */
unsigned long long power(unsigned long long x,
                                  long long y, long long p)
{
    unsigned long long res = 1; // Initialize result
  
    while (y > 0)
    {
     
        // If y is odd, multiply x with result
        if (y & 1)
            res = (res * x) % p;
 
        // y must be even now
        y = y >> 1; // y = y/2
        x = (x * x) % p;
    }
    return res;
}
 
// Returns n^(-1) mod p
unsigned long long modInverse(unsigned long long n, 
                                            long long p)
{
    return power(n, p - 2, p);
}
 
// Returns nCr % p using Fermat's little
// theorem.
unsigned long long nCrModPFermat(unsigned long long n,
                                 long long r, long long p)
{
    // If n<r, then nCr should return 0
    // Base case
    if (r == 0)
        return 1;
 
    // Fill factorial array so that we
    // can find all factorial of r, n
    // and n-r
    unsigned long long fac[n + 1];
    fac[0] = 1;
    for (int i = 1; i <= n; i++)
        fac[i] = (fac[i - 1] * i) % p;
 
    long long ncr = (fac[n] * modInverse(fac[r], p) % p
            * modInverse(fac[n - r], p) % p)
           % p;
    return ncr;
}
 
// Driver program
int main()
{
    // p must be a prime greater than n.
    int n = 41, r = 27, p = 139;
    cout << "Value of nCr % p is "
         << nCrModPFermat(n, r, p);
    return 0;
}
Python
def ncr(n, r, p):
    # initialize numerator
    # and denominator
    num = den = 1
    for i in range(r):
        num = (num * (n - i)) % p
        den = (den * (i + 1)) % p
    return (num * pow(den,
            p - 2, p)) % p
 
# p must be a prime
# greater than n
n, r, p = 41, 27, 139
print("Value of nCr % p is",
               ncr(n, r, p))

Output

Value of 41C27 % 139 is 78

Time Complexity:  O(n + log n) -> O(n).

Precomputing the factorial numbers from 1 to n will allow queries to be answered in O(n). And solving the power function will take Olog n) time complexity.

Space Complexity – O(n)

We are using an extra array to store all the values of factorial that we have precomputed  to solve nCr

Share your love

You may also like