Miller Rabin Probabilistic Primality Algorithm

Miller Rabin is undoubtedly the most popular probabilistic primality testing algorithm.

It uses the concept of Fermat’s Little Theorem.So if you don’t know fermat’s little theorem,it is advised to understand it first – Fermat

Let us Begin.

For an odd prime number p, we can write \mathbf{\left ( p-1 \right )= {2^{d}.s}}

Here s is odd and d>=1

By Fermat’s Little Theorem,we can write –

Thus for p to be prime , at least one of the terms in the expansion should be divisible by p.

In other words,

   or     or        or      ….. or .

If none of the conditions are true then is definitely composite, otherwise p maybe prime.

To increase the accuracy of the algorithm, we will perform the check for certain number of iterations for some random values of a lesser than p.

Here’s  the python implementation –

import random

def exp(a,p,m):
	if (p==0):
		return 1
	elif (p%2==1):
		return (a*exp(a,p-1,m))%m
	x=exp(a,p/2,m)
	return (x*x)%m

def miller_rabin(p):

	if (p==2 or p==3):return 1
	if (p%2==0):return 0
	if (p<3):return 0

	for i in range(10):
		#random 'a'
		a = random.randrange(2,p-2)

		s=p-1
		while(s%2==0):s/=2

		mod = exp(a,s,p)
		if (mod==1 or mod==p-1):
			continue
		
		flag=0
		while(s!=(p-1)):
			mod = exp(a,s,p)
			if (mod==(p-1)):
				flag=1
				break
			s*=2
		if (flag==0):
			return 0
	return 1

 

The Algorithm is slower than Fermat’s Little primality test but is still preferred due to it’s high accuracy.

NOTE – If you are implementing the algorithm in C/C++ , take care of the overflows while multiplying two long long ints  in modular exponentiation.

Try out this problem on SPOJ/PON

For any queries,drop a comment below 🙂

Advertisements

Fermat’s Little Theorem Proof & Algorithm

Fermat’s Little Theorem states that for a prime number p , the following statement holds true –

ap ≡ a (mod p)
(~ which can be read as ap mod p = a)
OR
ap-1 ≡ 1 (mod p)

Here is any number less than p.

PROOF

There are tons of proofs for Fermat’s Little Theorem on Wikipedia . I am discussing an easy one based on binomial expansion – Proof by Induction.

Consider a binomial expression – (1+a)p

On expanding it , we get –
(1+a)p = ap + pC1 ap-1 + pC2 ap-2 + pC3 ap-3 …… + 1

If we take modulo p both sides, All the terms except the first and the last term becomes 0.

How? Let’s See.

Consider a binomial coefficient – pC = p!(p-k)! (k)!

Since is prime (no divisors) ,pCk   will have p as it’s divisor .Thus modulo p will yield 0 as result.

we get –
(1+a)p ≡ ap + 1 mod p

Considering the Fermat’s Little theorem true , we can replace ap mod p by
We get the following result – (1+a)p ≡ (1+a) mod p

Thus Fermat’s Little Theorem is proved by Induction.

Now that we’ve proved the theorem , let’s see how we can use it to check the probabilistic primality of a number.

PRIME Algorithm

To check whether a number p is prime or not, just take a random a less than p and verify the theorem.

if  (ap-1 mod p)  yields 1 as result , then the number may be prime.

To increase the accuracy of the algorithm,perform the check for a certain number of iterations.

Here’s the python Code for the Algorithm

import random

def exp(a,p,m):
	if (p==0):
		return 1
	elif (p%2==1):
		return (a*exp(a,p-1,m))%m
	x=exp(a,p/2,m)
	return (x*x)%m

def fermat_primality(p):
	if (p==2 or p==3):return 1 
	if (p<=3):return 0

	#number of iterations = 10

	for i in range(10):	
		#Generating a random number a lesser than p
		a = random.randrange(2,p-1)	
		if (exp(a,p-1,p)!=1):
			return 0

	return 1

 

The time complexity of the Algorithm is O(iterations*log(n)) which is pretty fast. But the algorithm is not always accurate as there may be many composite numbers p which may yield the result for some value of a. Therefore it is advised to perform the check for certain number of iterations.

Also there are certain composite numbers n which satisfies the Fermat’s Little congruence relation for all values 1<a<n (a is relatively prime to n). Such numbers are called Carmichael numbers.

But they are very less in number – around 2000 under 25,000,0000,000.

Here’s a problem to try – SPOJ/PON

INVERSE MODULO

The Fermat’s result proved above ap-1 ≡ 1 (mod p)  can be rewritten as  a.ap-2 ≡ 1 (mod p).

In the expression a.x ≡ 1 (mod p) , x is the inverse of a.

Thus, when p is prime , (1a) mod p , is given by ap-2 mod p .

def exp(a,p,m):
	if (p==0):
		return 1
	elif (p%2==1):
		return (a*exp(a,p-1,m))%m
	x=exp(a,p/2,m)
	return (x*x)%m

def mod_inverse(a,p):
	return exp(a,p-2,p)

 

For any queries related to the article , feel free to drop a comment below 🙂

EDIT – If you’re trying to implement the above algorithm on C/C++, take care of overflows while multiplying two long long int numbers.