# 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 –

$\mathbf{a^{\left&space;(&space;p-1&space;\right&space;)}-1&space;\equiv&space;0&space;\pmod&space;p}$

$\mathbf{a^{2^{d}.s}-1&space;\equiv&space;0&space;\pmod&space;p}$

$\Rightarrow&space;\left&space;(&space;{a^{2^{d-1}.s}}&space;\right&space;)^2-1\equiv&space;0\pmod&space;p$

$\Rightarrow&space;\left&space;(&space;{a^{2^{d-1}.s}}&space;-1&space;\right&space;).\left&space;(&space;{a^{2^{d-1}.s}}&space;+1&space;\right&space;)\equiv&space;0\pmod&space;p$

$\Rightarrow&space;\left&space;(&space;{a^{2^{d-2}.s}}&space;-1&space;\right&space;).\left&space;(&space;{a^{2^{d-2}.s}}&space;+1&space;\right&space;).\left&space;(&space;{a^{2^{d-1}.s}}&space;+1&space;\right&space;)\equiv&space;0\pmod&space;p$

$\Rightarrow&space;\left&space;(&space;{a^{s}}&space;-1&space;\right&space;).\left&space;(&space;{a^{s}}&space;+1&space;\right&space;).\left&space;(&space;{a^{2s}}&space;+1&space;\right&space;).\left&space;(&space;{a^{4s}}&space;+1&space;\right&space;).......\left&space;(&space;{a^{2^{d-1}.s}}&space;+1&space;\right&space;)\equiv&space;0\pmod&space;p$

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

In other words,

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 🙂

# 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.

# Optimizing Sieve of Eratosthenes

One of the most popular prime number algorithm is the Sieve which comes handy when you are required to generate prime numbers up to an upper bound.In this post we’ll be optimizing sieve to save both memory as well as time.

The optimization that we’ll be performing will be both memory as well as time efficient.

NOTE – I won’t be discussing how the basic sieve works.(Check out – Sieve(without optimisation)).

We will move step by step towards optimizing the algorithm to maximum extent.(If you know any further optimizations, please do share in the comment box).

So here’s what sieve looks like without any optimization.

//============================================================================
// Name        : sieve.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

#define MAX 1000000
#define SQ 1000

int prime[MAX+1];

int sieve(){
for (int i=2;i<=SQ;i++){
if (!prime[i]){
for(int j=i*i;j<=MAX;j+=i){
prime[j] = 1;
}
}
}
return 0;
}

int main(){
sieve();
return 0;
}


Lets begin our optimization!

Considering that you are processing the algorithm for 1000 elements, the size of the array is 4000 bytes (4 bytes or 32 bits per int value).

Instead of using an entire integer(32 bits) for marking as prime(or composite) , you could use just 1 bit.How?

use array index ‘0’ to store values for the first 32 (from 1 to 32) numbers, index ‘1’ for the next 32 (from 33 to 64) numbers and so on..

This will save us a lot of memory.But still there’s no improvement in time.

Just don’t consider the even values at all.This will make your algorithm faster by a factor of 2 or more.

Now all we have are odd values to deal with.Lets see what the arrangement will look like after neglecting the even values.

This not only speeds up your algorithm, but saves even more memory.(Now we are using 4000/64 = 62.5 bytes of memory).

So here’s the code for the complete optimization and later we’ll see how the bit wise magic is working.

//============================================================================
// Name        : optimised_bitwise_sieve.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

#define MAX 100000000
#define SQ 10000

//These macros are discussed below
#define check(n) (prime[n>>6]&(1<<((n&63)>>1)))
#define set(n) prime[n>>6]|=(1<<((n&63)>>1))

int prime[MAX>>6];

int sieve(){
for(int i=3;i<=SQ;i+=2){
if (!check(i)){
int inc = 2*i;
for(int j=i*i;j<=MAX;j+=inc){
set(j);
}
}
}
return 0;
}

int main(){
sieve();
return 0;
}


Let us see how we can access and modify the bit value for a number ‘n’ (n is odd) .

1. Get the array index where the bit value of ‘n’ is stored. Since we are storing the bit values of 32 consecutive odd numbers in one integer , the array index can be generalized as n/64 or n>>6 . So the array value if prime[n>>6].
2. Now in order see which bit to visit,we have to modulo by 64 and then divide by 2.Example – numbers 1 and 65 in the picture above both have their result stored in first (least significant) bit but in different indexes. If we modulo both of them by 64,we will get 1 and then divide by , we will get which is the required bit position.
3. the modulo operation n%64 can be seen as n&63 . This is similar to what happens when we divide a number(in base 10) by 10^k . We get the last k digits of that number.Except this time the modulo operation is taking place in base 2. So if we modulo n number by 2^6 , we will get the last 6 digits in the base 2 representation of n. (which is n&63).Try it out yourself for further clarification.

So to sum up, the first point gives us the logic for the array index in which the bit value for n is present – prime[n>>6]

And second and third points combined tells us about the bit position given by – ((n&63)>>1)

Now we can check as well as set the bit values.

Hope you understood the post 🙂

Here’s a question on optimized bit wise sieve (where normal sieve won’t work) on SPOJ-TDKPRIME.

# Tweaking sort() function C++

C++ comes with inbuilt sort() function in the ‘algorithm.h’ module which can be really handy to both beginner as well as experienced programmers.In this post,we’ll try to explore sort() function and extend it’s functionality by altering the comparison function in the sorting algorithm.

Lets get started!

Syntax – sort(first,last,comp)

• Here,’first’ and ‘last’ refers to the pointer to the first and last element.
• ‘comp’ is the name of the function that accepts two values,compares them and then returns a boolean value.it is absolutely optional.
• In the absence of ‘comp’ , an operator overloaded ‘<‘ function name is used.

So,lets discuss two techniques to use sort function.

//============================================================================
// Name        : sort_func.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

struct node{
int val;
};

node a[5];

//this function only accepts struct/class type arguments.

bool operator <(node x,node y){
//if x.val is less than y.val return true else false
return x.val<y.val;
}

int main(){
a[0].val = 2;
a[1].val = 1;
a[2].val = 4;
a[3].val = 3;
a[4].val = 5;

sort(a,a+5);
for(int i=0;i<5;i++) cout << a[i].val << " ";

return 0;
}


Output – 1 2 3 4 5

Another technique is by providing a name to the function.We’ll take the name ‘func’.

//============================================================================
// Name        : sort_func.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

vector<int> a;

bool func(int x,int y){
return x<y;
}

int main(){
a.push_back(2);
a.push_back(1);
a.push_back(4);
a.push_back(3);
a.push_back(5);

sort(a.begin(),a.end(),func);

for(int i=0;i<5;i++) cout << a[i] << " ";

return 0;
}


Output – 1 2 3 4 5

Since you are using a custom function , you may pass arguments of any but same type.Also this time it is mandatory to pass the function name as the third argument.

Note,here the function is able to accept ‘int’ values(since we are using a custom function name).If we use the operator overloaded function(first technique) , then we are bound to use struct/class only.

Now , lets see how the sort() function behaves if we alter the comparison function.

Lets begin by changing the direction of the comparison operator.

//============================================================================
// Name        : sort_func.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

bool comp(int x,int y){
return x>y;	//returns true if x>y(reversed the order)
}

int main(){
int x[] = {2,1,4,3,5};
vector<int> a(x,x+5);
sort(a.begin(),a.end(),comp);
for(int i=0;i<5;i++) cout << a[i] << " ";
return 0;
}

//OUTPUT - 5 4 3 2 1


This time we’ve got the result printed in non-increasing format.

And here’s how you can sort strings .

//============================================================================
// Name        : sort_func.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

struct node{
char s[10];
};

node a[5];

bool operator <(node x,node y){
if (strcmp(x.s,y.s)<0) return true;
else return false;
}

int main(){
strcpy(a[0].s,"aaca");
strcpy(a[1].s,"aca");
strcpy(a[2].s,"aaaa");
strcpy(a[4].s,"baca");
sort(a,a+5);
for(int i=0;i<5;i++) cout << a[i].s << " ";
return 0;
}

//OUTPUT - aaaa aaca aca baca zadd


Now suppose, You have to sort various cartesian co-ordinates given such that smaller x co-ordinate should appear before larger x co-ordinate.If two co-ordinates have same value of x , then sort according to y.

Here’s the code –

//============================================================================
// Name        : sort_func.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <bits/stdc++.h>
using namespace std;

struct node{
int x,y;
};

node c[5];

bool operator <(node a,node b){
if (a.x==b.x) return a.y<b.y;
else return a.x<b.x;
}

int main(){
c[0].x = 1; c[0].y = 2;
c[1].x = 1; c[1].y = 4;
c[2].x = 3; c[2].y = 7;
c[3].x = 1; c[3].y = 1;
c[4].x = 2; c[4].y = 7;
sort(c,c+5);
for(int i=0;i<5;i++) cout << "(" << c[i].x << "," << c[i].y << ")" << " ";
return 0;
}

//INPUT  - (1,2) (1,4) (3,7) (1,1) (2,7)
//OUTPUT - (1,1) (1,2) (1,4) (2,7) (3,7)


Depending upon different situations, modifying the comparison function can help in making your code compact as well as speeding up your programming.

Feel free to drop a comment below 🙂

# Mark next greater element in an array.

This is an interesting problem I found on geeksforgeeks .

### Problem Statement

You are given an array of integers. For each element in that array, you need to find the position of next element greater than it.

### Analysis

The above problem can easily be solved in O(n*n) using the naive algorithm. The naive algorithm involves checking the greater element for each element in O(n). Thus , the complexity for the entire algorithm becomes O(n*n).

//============================================================================
// Name        : Greater_element_naive.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <iostream>
#include <string.h>
using namespace std;

#define MAX 5000

int n,a[MAX];

int solve(){
for (int i=0;i<n;i++){
int flag = 0;	// flag to check if there exist an element greater and to the right of i
for(int j=i+1;j<n;j++){
if (a[j]>a[i]){
flag=1;
cout << i << ":" << j << endl;
break;
}
}
if (flag==0) cout << i << ":" << -1 << endl;
}
return 0;
}

int input(){
cin >> n;
for(int i=0;i<n;i++) cin >> a[i];
return 0;
}

int main() {
input();
solve();
return 0;
}


However the above algorithm is too slow and requires optimization.

So,Here’s a O(n) approach using stacks.

### Procedure

Start by traversing the array and maintaining a stack. Initially the stack is empty denoted by top = -1.

You will push an element and increment the array pointer into the stack if –

• The stack is empty
• The element in the array is lesser than the top element of the stack.

You will pop the element (but do not increment the array pointer) if –

•   The element in the array is greater than the top element of the stack.(This time you’ve found the element greater for the top element of the stack therefore you also mark the result).

//============================================================================
// Name        : Greater_element_optimised.cpp
// Author      : Prateek Agarwal
// Institution : NIT Jalandhar
//============================================================================

#include <iostream>
#include <string.h>
using namespace std;

#define MAX 5000

int n,a[MAX];
int mark[MAX];

struct stack{
int value;
int index; //for storing the index of the element in array 'a' pushed into the stack
};
stack s[MAX];
int top = -1;

//O(n) approach using stack

int solve(){
int i=0;
while(i<n){
if (top==-1 || s[top].value >= a[i]){
top++;
s[top].value = a[i];	//pushing the element into the stack
s[top].index = i;
i++;
}
else{
mark[s[top].index] = i;
top--;	//deleting an element from the top
}
}

//displaying the result

for (int i=0;i<n;i++){
cout << i << ":" << mark[i] << endl;
}

return 0;
}

int input(){
memset(mark,-1,sizeof(mark));
cin >> n;
for(int i=0;i<n;i++) cin >> a[i];
return 0;
}

int main() {
input();
solve();
return 0;
}


### Example

Consider the following input – n = 4 , a = [3,2,4,1]

1. Initially the stack is empty. So , top = -1.Also i(array pointer) = 0.
2. Since the stack is empty,push the first element and increment the i. The top element in the stack becomes 3.
3. The next element i.e 2 is lesser than the top element of the stack i.e 3. So,push it as well and increment i.The top element in the stack becomes 2.
4. The next element i.e 4 is greater than the top element of the stack. We’ve got the greater element for top element of the stack.Therefore mark the result and pop out the top element.(Do not increment array pointer).The top element is now 3.
5. the array element 4 is greater than the top element of the stack i.e 3.Therefore mark the result pop the top element and continue.
6. Since the stack is empty now,push the array element 4.increment array pointer.
7. the next element is 1 (smaller than top element of stack).So,push it and increment i.
8. i=n , End of loop.
9. Since the last two elements in the array aren’t marked in the loop,it means that there exists no such result for these values.(The default  values are set as -1).

The above code produces the following output for the sample example-

0:2
1:2
2:-1
3:-1

### Analysis of Complexity – O(n)

Consider any case, the above code involves pushing and popping each element into the stack exactly once.

Thus, the complexity is O(n) which is a big improvement over the previously discussed naive algorithm.

For any query , drop a comment below.