Open In App

Compute nCr%p using Fermat Little Theorem

Improve
Improve
Like Article
Like
Save
Share
Report

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

Input:  n = 10, r = 2, p = 13
Output: 6
Explanation: 10C2 is 45 and 45 % 13 is 6.

Input:  n = 6, r = 2, p = 13
Output: 2
Recommended Practice

We have discussed the following methods in previous posts. 
Compute nCr % p | Set 1 (Introduction and Dynamic Programming Solution) 
Compute nCr % p | Set 2 (Lucas Theorem)
In this post, Fermat Theorem-based solution is discussed.
Background: 
Fermat’s little theorem and modular inverse 
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) 
For example, if a = 2 and p = 7, 27 = 128, and 128 – 2 = 7 × 18 is an integer multiple of 7.
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, i.e 
ap-1 = 1 (mod p)
If we multiply both sides by a-1, we get. 
ap-2 = a-1 (mod p)
So we can find modular inverse as p-2
Computation: 

We know the formula for  nCr 
nCr = fact(n) / (fact(r) x fact(n-r)) 
Here fact() means factorial.

 nCr % p = (fac[n]* modIverse(fac[r]) % p *
               modIverse(fac[n-r]) % p) % p;
Here modIverse() means modular inverse under
modulo p.

Following is the implementation of the above algorithm. In the following implementation, an array fac[] is used to store all the computed factorial values. 

C++




// A modular inverse based solution to
// compute nCr % p
#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,
                                  int y, int p)
{
    unsigned long long res = 1; // Initialize result
 
    x = x % p; // Update x if it is more than or
    // equal to p
 
    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, 
                                            int p)
{
    return power(n, p - 2, p);
}
 
// Returns nCr % p using Fermat's little
// theorem.
unsigned long long nCrModPFermat(unsigned long long n,
                                 int r, int p)
{
    // If n<r, then nCr should return 0
    if (n < r)
        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;
 
    return (fac[n] * modInverse(fac[r], p) % p
            * modInverse(fac[n - r], p) % p)
           % p;
}
 
// Driver program
int main()
{
    // p must be a prime greater than n.
    int n = 10, r = 2, p = 13;
    cout << "Value of nCr % p is "
         << nCrModPFermat(n, r, p);
    return 0;
}


Java




// A modular inverse based solution to
// compute nCr %
import java.io.*;
 
class GFG {
 
    /* Iterative Function to calculate
    (x^y)%p in O(log y) */
    static int power(int x, int y, int p)
    {
 
        // Initialize result
        int 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;
    }
 
    // Returns n^(-1) mod p
    static int modInverse(int n, int p)
    {
        return power(n, p - 2, p);
    }
 
    // Returns nCr % p using Fermat's
    // little theorem.
    static int nCrModPFermat(int n, int r,
                             int p)
    {
 
          if (n<r)
              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
        int[] fac = new int[n + 1];
        fac[0] = 1;
 
        for (int i = 1; i <= n; i++)
            fac[i] = fac[i - 1] * i % p;
 
        return (fac[n] * modInverse(fac[r], p)
                % p * modInverse(fac[n - r], p)
                % p)
            % p;
    }
 
    // Driver program
    public static void main(String[] args)
    {
 
        // p must be a prime greater than n.
        int n = 10, r = 2, p = 13;
        System.out.println("Value of nCr % p is "
                           + nCrModPFermat(n, r, p));
    }
}
 
// This code is contributed by Anuj_67.


Python3




# Python3 program to calculate nCr % p
 
#Python function to calculate nCr % p
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 = 10, 2, 13
print("Value of nCr % p is", ncr(n, r, p))


C#




// A modular inverse based solution to
// compute nCr % p
using System;
 
class GFG {
 
    /* Iterative Function to calculate
    (x^y)%p in O(log y) */
    static int power(int x, int y, int p)
    {
 
        // Initialize result
        int 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;
    }
 
    // Returns n^(-1) mod p
    static int modInverse(int n, int p)
    {
        return power(n, p - 2, p);
    }
 
    // Returns nCr % p using Fermat's
    // little theorem.
    static int nCrModPFermat(int n, int r,
                             int p)
    {
 
      if (n<r)
            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
        int[] fac = new int[n + 1];
        fac[0] = 1;
 
        for (int i = 1; i <= n; i++)
            fac[i] = fac[i - 1] * i % p;
 
        return (fac[n] * modInverse(fac[r], p)
                % p * modInverse(fac[n - r], p)
                % p)
            % p;
    }
 
    // Driver program
    static void Main()
    {
 
        // p must be a prime greater than n.
        int n = 10, r = 11, p = 13;
        Console.Write("Value of nCr % p is "
                      + nCrModPFermat(n, r, p));
    }
}
 
// This code is contributed by Anuj_67


PHP




<?php
// A modular inverse
// based solution to
// compute nCr % p
 
// Iterative Function to
// calculate (x^y)%p
// in O(log y)
function power($x, $y, $p)
{
     
    // Initialize result
    $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 & 1)
            $res = ($res * $x) % $p;
 
        // y must be
        // even now
        // y = y/2
        $y = $y >> 1;
        $x = ($x * $x) % $p;
    }
    return $res;
}
 
// Returns n^(-1) mod p
function modInverse($n, $p)
{
    return power($n, $p - 2, $p);
}
 
// Returns nCr % p using
// Fermat's little
// theorem.
function nCrModPFermat($n, $r, $p)
{
     
    if ($n<$r)
          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
    //$fac[$n+1];
    $fac[0] = 1;
    for ($i = 1; $i <= $n; $i++)
        $fac[$i] = $fac[$i - 1] *
                        $i % $p;
 
    return ($fac[$n] * modInverse($fac[$r], $p) % $p *
             modInverse($fac[$n - $r], $p) % $p) % $p;
}
 
    // Driver Code
    // p must be a prime
    // greater than n.
    $n = 10;
    $r = 2;
    $p = 13;
    echo "Value of nCr % p is ",
         nCrModPFermat($n, $r, $p);
         
// This code is contributed by Ajit.
?>


Javascript




<script>
 
    // A modular inverse based solution
    // to compute nCr % p
     
    /* Iterative Function to calculate
    (x^y)%p in O(log y) */
    function power(x, y, p)
    {
  
        // Initialize result
        let 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;
    }
  
    // Returns n^(-1) mod p
    function modInverse(n, p)
    {
        return power(n, p - 2, p);
    }
  
    // Returns nCr % p using Fermat's
    // little theorem.
    function nCrModPFermat(n, r, p)
    {
  
      if (n<r)
      {
          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
      let fac = new Array(n + 1);
      fac[0] = 1;
 
      for (let i = 1; i <= n; i++)
        fac[i] = fac[i - 1] * i % p;
 
      return (fac[n] * modInverse(fac[r], p) % p *
              modInverse(fac[n - r], p) % p) % p;
    }
     
    // p must be a prime greater than n.
    let n = 10, r = 2, p = 13;
    document.write("Value of nCr % p is " +
                    nCrModPFermat(n, r, p));
 
</script>


Output

Value of nCr % p is 6

Time Complexity: O(n)
Auxiliary Space: O(n)

We can further improve it’s space complexity:

Instead of calculating factorial in a different array, we can directly multiply numbers with some cancellations.

We know the formula for  nCr 
nCr = fact(n) / (fact(r) x fact(n-r)) 
fact(n) = n * (n-1) * (n-2) * (n-3)* .... * 1;
fact(r) = r * (r-1) * (r-2) * ......... *1;
Because r is always less than n in nCr So all factors in fact(r) also come in fact(n),
hence we can cancell them.
And get the multiplication of rest elements of numerator and denominator(fact(n-r))
Note that rest elements of numerator will be n* (n-1) *(n-2) * .... (r+1).

We can also improve time complexity with the logic nCr is equal to nC(n-r). So whenever n-r is less than r compute nCn-r instead of nCr.

See the below implementation:

C++




// A modular inverse based solution to
// compute nCr % p
#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, int y, int p)
{
    unsigned long long res = 1; // Initialize result
 
    x = x % p; // Update x if it is more than or
    // equal to p
 
    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, int p)
{
    return power(n, p - 2, p);
}
unsigned long long mul(unsigned long long x,
                       unsigned long long y, int p)
{
    return x * 1ull * y % p;
}
unsigned long long divide(unsigned long long x,
                          unsigned long long y, int p)
{
    return mul(x, modInverse(y, p), p);
}
// Returns nCr % p using Fermat's little
// theorem.
unsigned long long nCrModPFermat(unsigned long long n,
                                 int r, int p)
{
    // If n<r, then nCr should return 0
    if (n < r)
        return 0;
    // Base case
    if (r == 0)
        return 1;
    // if n-r is less calculate nCn-r
    if (n - r < r)
        return nCrModPFermat(n, n - r, p);
 
    // Fill factorial array so that we
    // can find all factorial of r, n
    // and n-r
    unsigned long long res = 1;
    // keep multiplying numerator terms and dividing denominator terms in res
    for (int i = r; i >= 1; i--)
        res = divide(mul(res, n - i + 1, p), i, p);
    return res;
}
 
// Driver program
int main()
{
    // p must be a prime greater than n.
    int n = 10, r = 2, p = 13;
    cout << "Value of nCr % p is "
         << nCrModPFermat(n, r, p);
    return 0;
}
 
//Code and idea by Harsh Singh (hsnooob)


Java




import java.util.*;
 
public class Gfg {
 
  // Iterative Function to calculate (x^y)%p in O(log y)
  public static long power(long x, int y, int p) {
    long res = 1; // Initialize result
 
    x = x % p; // Update x if it is more than or equal to p
 
    while (y > 0) {
      // If y is odd, multiply x with result
      if ((y & 1) == 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
  public static long modInverse(long n, int p) {
    return power(n, p - 2, p);
  }
 
  public static long mul(long x, long y, int p) {
    return x * 1L * y % p;
  }
 
  public static long divide(long x, long y, int p) {
    return mul(x, modInverse(y, p), p);
  }
 
  // Returns nCr % p using Fermat's little theorem.
  public static long nCrModPFermat(long n, long r, int p) {
    // If n<r, then nCr should return 0
    if (n < r)
      return 0;
 
    // Base case
    if (r == 0)
      return 1;
 
    // if n-r is less calculate nCn-r
    if (n - r < r)
      return nCrModPFermat(n, n - r, p);
 
    // Fill factorial array so that we can find all factorial of r, n and n-r
    long res = 1;
    // keep multiplying numerator terms and dividing denominator terms in res
    for (long i = r; i >= 1; i--)
      res = divide(mul(res, n - i + 1, p), i, p);
    return res;
  }
 
  // Driver program
  public static void main(String[] args) {
    // p must be a prime greater than n.
    int n = 10, r = 2, p = 13;
    System.out.println("Value of nCr % p is " + nCrModPFermat(n, r, p));
  }
}


Python3




def power(x, y, p):
    res = 1
    x = x % p
 
    while y > 0:
        if y & 1:
            res = (res * x) % p
        y = y >> 1
        x = (x * x) % p
    return res
 
def modInverse(n, p):
    return power(n, p - 2, p)
 
def mul(x, y, p):
    return (x * y) % p
 
def divide(x, y, p):
    return mul(x, modInverse(y, p), p)
 
def nCrModPFermat(n, r, p):
    if n < r:
        return 0
    if r == 0:
        return 1
    if n - r < r:
        return nCrModPFermat(n, n - r, p)
    res = 1
    for i in range(1, r + 1):
        res = divide(mul(res, n - i + 1, p), i, p)
    return res
 
n = 10
r = 2
p = 13
print("Value of nCr % p is", nCrModPFermat(n, r, p))


C#




using System;
 
class Program {
    // Iterative Function to calculate (x^y)%p in O(log y)
    static long Power(long x, int y, int p)
    {
        long res = 1; // Initialize result
 
        x = x % p; // Update x if it is more than or equal
                   // to p
 
        while (y > 0) {
            // If y is odd, multiply x with result
            if ((y & 1) == 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
    static long ModInverse(long n, int p)
    {
        return Power(n, p - 2, p);
    }
 
    static long Mul(long x, long y, int p)
    {
        return x * 1L * y % p;
    }
 
    static long Divide(long x, long y, int p)
    {
        return Mul(x, ModInverse(y, p), p);
    }
 
    // Returns nCr % p using Fermat's little theorem.
    static long NCrModPFermat(long n, long r, int p)
    {
        // If n<r, then nCr should return 0
        if (n < r)
            return 0;
 
        // Base case
        if (r == 0)
            return 1;
 
        // if n-r is less calculate nCn-r
        if (n - r < r)
            return NCrModPFermat(n, n - r, p);
 
        // Fill factorial array so that we can find all
        // factorial of r, n and n-r
        long res = 1;
        // keep multiplying numerator terms and dividing
        // denominator terms in res
        for (long i = r; i >= 1; i--)
            res = Divide(Mul(res, n - i + 1, p), i, p);
        return res;
    }
 
    static void Main(string[] args)
    {
        // p must be a prime greater than n.
        int n = 10, r = 2, p = 13;
        Console.WriteLine("Value of nCr % p is "
                          + NCrModPFermat(n, r, p));
    }
}


Javascript




// A modular inverse based solution to
// compute nCr % p
 
function power(x, y, p) {
  let res = 1; // Initialize result
 
  x = x % p; // Update x if it is more than or equal to p
 
  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
function modInverse(n, p) {
  return power(n, p - 2, p);
}
 
function mul(x, y, p) {
  return (x * y) % p;
}
 
function divide(x, y, p) {
  return mul(x, modInverse(y, p), p);
}
 
// Returns nCr % p using Fermat's little theorem.
function nCrModPFermat(n, r, p) {
  // If n<r, then nCr should return 0
  if (n < r) return 0;
  // Base case
  if (r == 0) return 1;
  // if n-r is less calculate nCn-r
  if (n - r < r) return nCrModPFermat(n, n - r, p);
 
  // Fill factorial array so that we
  // can find all factorial of r, n
  // and n-r
  let res = 1;
  // keep multiplying numerator terms and dividing denominator terms in res
  for (let i = r; i >= 1; i--)
    res = divide(mul(res, n - i + 1, p), i, p);
  return res;
}
 
// Driver program
let n = 10,
  r = 2,
  p = 13;
console.log("Value of nCr % p is " + nCrModPFermat(n, r, p));


Output

Value of nCr % p is 6

Time Complexity: O(n)
Auxiliary Space: O(1)

Improvements: 
In competitive programming, we can pre-compute fac[] for a given upper limit so that we don’t have to compute it for every test case. We also can use unsigned long long int everywhere to avoid overflows.



Last Updated : 15 Apr, 2023
Like Article
Save Article
Previous
Next
Share your thoughts in the comments
Similar Reads