A blog about programming and more programming.

Posts Tagged ‘fermat

Modular Multiplicative Inverse

with 35 comments

The modular multiplicative inverse of an integer a modulo m is an integer x such that a^{-1} \equiv x \pmod{m}.

That is, it is the multiplicative inverse in the ring of integers modulo m. This is equivalent to ax \equiv aa^{-1} \equiv 1 \pmod{m}.

The multiplicative inverse of a modulo m exists if and only if a and m are coprime (i.e., if gcd(a, m) = 1).

Let’s see various ways to calculate Modular Multiplicative Inverse:

1. Brute Force
We can calculate the inverse using a brute force approach where we multiply a with all possible values x and find a x such that ax \equiv 1 \pmod{m}. Here’s a sample C++ code:

int modInverse(int a, int m) {
	a %= m;
	for(int x = 1; x < m; x++) {
		if((a*x) % m == 1) return x;

The time complexity of the above codes is O(m).

2. Using Extended Euclidean Algorithm
We have to find a number x such that a·x = 1 (mod m). This can be written as well as a·x = 1 + m·y, which rearranges into a·x – m·y = 1. Since x and y need not be positive, we can write it as well in the standard form, a·x + m·y = 1.

In number theory, Bézout’s identity for two integers a, b is an expression ax + by = d, where x and y are integers (called Bézout coefficients for (a,b)), such that d is a common divisor of a and b. If d is the greatest common divisor of a and b then Bézout’s identity ax + by = gcd(a,b) can be solved using Extended Euclidean Algorithm.

The Extended Euclidean Algorithm is an extension to the Euclidean algorithm. Besides finding the greatest common divisor of integers a and b, as the Euclidean algorithm does, it also finds integers x and y (one of which is typically negative) that satisfy Bézout’s identity
ax + by = gcd(a,b). The Extended Euclidean Algorithm is particularly useful when a and b are coprime, since x is the multiplicative inverse of a modulo b, and y is the multiplicative inverse of b modulo a.

We will look at two ways to find the result of Extended Euclidean Algorithm.

Iterative Method
This method computes expressions of the form ri = axi + byi for the remainder in each step i of the Euclidean algorithm. Each successive number ri can be written as the remainder of the division of the previous two such numbers, which remainder can be expressed using the whole quotient qi of that division as follows:
r_i = r_{i-2} - q_i r_{i-1}.\,
By substitution, this gives:
r_i = (ax_{i-2} + by_{i-2}) - q_i (ax_{i-1} + by_{i-1}),\,which can be written
r_i = a(x_{i-2} - q_i x_{i-1}) + b(y_{i-2} - q_i y_{i-1}).\,
The first two values are the initial arguments to the algorithm:
r_1 = a = a\times1 + b\times0\,
r_2 = b = a\times0 + b\times1.\,
So the coefficients start out as x1 = 1, y1 = 0, x2 = 0, and y2 = 1, and the others are given by
x_i = x_{i-2} - q_i x_{i-1},\,
y_i = y_{i-2} - q_i y_{i-1}.\,
The expression for the last non-zero remainder gives the desired results since this method computes every remainder in terms of a and b, as desired.

So the algorithm looks like,

  1. Apply Euclidean algorithm, and let qn(n starts from 1) be a finite list of quotients in the division.
  2. Initialize x0, x1 as 1, 0, and y0, y1 as 0,1 respectively.
    1. Then for each i so long as qi is defined,
    2. Compute xi+1 = xi-1qixi
    3. Compute yi+1 = yi-1qiyi
    4. Repeat the above after incrementing i by 1.
  3. The answers are the second-to-last of xn and yn.
/* This function return the gcd of a and b followed by
	the pair x and y of equation ax + by = gcd(a,b)*/
pair<int, pair<int, int> > extendedEuclid(int a, int b) {
	int x = 1, y = 0;
	int xLast = 0, yLast = 1;
	int q, r, m, n;
	while(a != 0) {
		q = b / a;
		r = b % a;
		m = xLast - q * x;
		n = yLast - q * y;
		xLast = x, yLast = y;
		x = m, y = n;
		b = a, a = r;
	return make_pair(b, make_pair(xLast, yLast));

int modInverse(int a, int m) {
    return (extendedEuclid(a,m).second.first + m) % m;

Recursive Method
This method attempts to solve the original equation directly, by reducing the dividend and divisor gradually, from the first line to the last line, which can then be substituted with trivial value and work backward to obtain the solution.
Notice that the equation remains unchanged after decomposing the original dividend in terms of the divisor plus a remainder, and then regrouping terms. So the algorithm looks like this:

  1. If b = 0, the algorithm ends, returning the solution x = 1, y = 0.
  2. Otherwise:
    • Determine the quotient q and remainder r of dividing a by b using the integer division algorithm.
    • Then recursively find coefficients s, t such that bs + rt divides both b and r.
    • Finally the algorithm returns the solution x = t, and y = s − qt.

Here’s a C++ implementation:

/* This function return the gcd of a and b followed by
	the pair x and y of equation ax + by = gcd(a,b)*/
pair<int, pair<int, int> > extendedEuclid(int a, int b) {
	if(a == 0) return make_pair(b, make_pair(0, 1));
	pair<int, pair<int, int> > p;
	p = extendedEuclid(b % a, a);
	return make_pair(p.first, make_pair(p.second.second - p.second.first*(b/a), p.second.first));

int modInverse(int a, int m) {
    return (extendedEuclid(a,m).second.first + m) % m;

The time complexity of the above codes is O(log(m)2).

3. Using Fermat’s Little Theorem
Fermat’s little theorem states that if m is a prime and a is an integer co-prime to m, then ap − 1 will be evenly divisible by m. That is a^{m-1} \equiv 1 \pmod{m}. or a^{m-2} \equiv a^{-1} \pmod{m}. Here’s a sample C++ code:

/* This function calculates (a^b)%MOD */
int pow(int a, int b, int MOD) {
int x = 1, y = a;
	while(b > 0) {
		if(b%2 == 1) {
			if(x>MOD) x%=MOD;
		y = (y*y);
		if(y>MOD) y%=MOD;
		b /= 2;
	return x;

int modInverse(int a, int m) {
    return pow(a,m-2,m);

The time complexity of the above codes is O(log(m)).

4. Using Euler’s Theorem
Fermat’s Little theorem can only be used if m is a prime. If m is not a prime we can use Euler’s Theorem, which is a generalization of Fermat’s Little theorem. According to Euler’s theorem, if a is coprime to m, that is, gcd(a, m) = 1, then a^{\varphi(m)} \equiv 1 \pmod{m}, where where φ(m) is Euler Totient Function. Therefore the modular multiplicative inverse can be found directly: a^{\varphi(m)-1} \equiv a^{-1} \pmod{m}. The problem here is finding φ(m). If we know φ(m), then it is very similar to above method.

Now lets take a little different question. Now suppose you have to calculate the inverse of first n numbers. From above the best we can do is O(n log(m)). Can we do any better? Yes.

We can use sieve to find a factor of composite numbers less than n. So for composite numbers inverse(i) = (inverse(i/factor(i)) * inverse(factor(i))) % m, and we can use either Extended Euclidean Algorithm or Fermat’s Theorem to find inverse for prime numbers. But we can still do better.

a * (m / a) + m % a = m
(a * (m / a) + m % a) mod m = m mod m, or
(a * (m / a) + m % a) mod m = 0, or
(- (m % a)) mod m = (a * (m / a)) mod m.
Dividing both sides by (a * (m % a)), we get
– inverse(a) mod m = ((m/a) * inverse(m % a)) mod m
inverse(a) mod m = (- (m/a) * inverse(m % a)) mod m

Here’s a sample C++ code:

vector<int> inverseArray(int n, int m) {
	vector<int> modInverse(n + 1,0);
	modInverse[1] = 1;
	for(int i = 2; i <= n; i++) {
		modInverse[i] = (-(m/i) * modInverse[m % i]) % m + m;
	return modInverse;

The time complexity of the above code is O(n).



Written by fR0DDY

October 9, 2011 at 12:29 AM

Miller Rabin Primality Test

with 17 comments

Miller Rabin Primality Test is a probabilistic test to check whether a number is a prime or not. It relies on an equality or set of equalities that hold true for prime values, then checks whether or not they hold for a number that we want to test for primality.


1> Fermat’s little theorem states that if p is a prime and 1 ≤ a < p then a^{p-1}\equiv1\pmod{p}.
2> If p is a prime and x^2\equiv1\pmod{p} or \left(x-1\right)\left(x+1\right)\equiv0\pmod{p} then x\equiv1\pmod{p} or x\equiv-1\pmod{p}.
3> If n is an odd prime then n-1 is an even number and can be written as 2^{s}.d. By Fermat’s Little Theorem either a^{d}\equiv1\pmod{n} or a^{2^r\cdot d}\equiv -1\pmod{n} for some 0 ≤ r ≤  s-1.
4> The Miller–Rabin primality test is based on the contrapositive of the above claim. That is, if we can find an a such that a^{d}\not\equiv1\pmod{n} and a^{2^r\cdot d}\not\equiv -1\pmod{n} for all 0 ≤ r ≤  s-1 then a is witness of compositeness of n and we can say n is not a prime. Otherwise, n may be a prime.
5> We test our number N for some random a and either declare that N is definitely a composite or probably a prime. The probably that a composite number is returned as prime after k itereations is 4^{-k}.


Input :A number N to be tested and a variable iteration-the number 
of 'a' for which algorithm will test N.
Output :0 if N is definitely a composite and 1 if N is probably a prime.

Write N as 2^{s}.d
For each iteration
	Pick a random a in [1,N-1]
	x = a^{d} mod n
	if x =1 or x = n-1 
		Next iteration
	for r = 1 to s-1
		x  = x^{2} mod n
		if x = 1
			return false
		if x = N-1
			Next iteration
	return false
return true

Here’s a python implementation :

import random
def modulo(a,b,c):
        x = 1
        y = a
        while b>0:
                if b%2==1:
                        x = (x*y)%c
                y = (y*y)%c
                b = b/2
        return x%c
def millerRabin(N,iteration):
        if N<2:
                return False
        if N!=2 and N%2==0:
                return False
        while d%2==0:
                d = d/2
        for i in range(iteration):
                a = random.randint(1, N-1)
                temp = d
                x = modulo(a,temp,N)
                while (temp!=N-1 and x!=1 and x!=N-1):
                        x = (x*x)%N
                        temp = temp*2
                if (x!=N-1 and temp%2==0):
                        return False
        return True


Written by fR0DDY

September 18, 2010 at 12:23 PM

%d bloggers like this: