Compute nCr%p using Fermat Little Theorem
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
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> |
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)); |
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.
Contact Us