# Checking if a number is prime or non prime

#### Miller-Rabin’s primality tests for checking if a number is prime or non prime.

Miller-Rabin’s method checks how likely it is that a given number is prime or non-prime. Thus, this primality test is probabilistic in nature.

Concepts to know before proceeding with Miller-Rabin’s primarlity test.

#### Fermat’s Theorem

If p is a prime number then for any integer a, the below Fermat’s test always hold true
a p - 1 = 1 ( mod p ) for all a ∈ { 1, 2, . . . , p − 1 } i.e it basically means that if we divide a p - 1 by p we would get a remainder of 1.
or
a p - 1 - 1 = 0 ( mod p )

Example 1: Say p = 5 and a = 2, then 2 5-1 = 2 4 = 16 % 5 = 1
Example 2: Say p = 7 and a = 3, then 3 7-1 = 3 6 = 729 % 7 = 1

Thus a prime number would always pass Fermat’s test.

• i.e If p is prime then a p - 1 - 1 = 0 ( mod p ) for all a ∈ { 1, 2, . . . , p − 1 }.

But the reverse is not true.

• i.e If a p - 1 - 1 = 0 ( mod p ) for all a ∈ { 1, 2, . . . , p − 1 } does not necessarily mean that p is prime.

Conclusion : A passed Fermat’s primality test would not tell us with certanity if a number is prime.

#### Miller-Rabin Primality Test

Some essential concepts

Equation for factorizing difference of squares is
x2 - y2 = ( x + y ) . ( x - y )

Powers

• ( a m ) n = ( a n ) m = a m . n
• a m . a n = a m + n
• a d . a d = a 2d
• a 2d . a 2d = a 4d = a 2 2 . d

Fermat’s equation for checking if a number p is prime is :
a p - 1 - 1 = 0 ( mod p ) for all a ∈ { 1, 2, . . . , p − 1 }.

As p is odd, we can write it as p - 1 = 2 r . d, ( d being odd ) so that we could factorize the equation using the difference of squares.
Thus we get the below equations

=> ( a 2 r - 1 . d - 1 ) . ( a 2 r - 1 . d + 1 ) = 0 ( mod p )

If 2 r - 1 is even, it could be written as 2 r - 1 = 2 . 2 r - 2 , thus we can further factorize the above equation as below else we stop factorizing

=> ( a 2 r - 2 . d - 1 ) . ( a 2 r - 2 . d + 1 ) . ( a 2 r - 1 . d + 1 ) = 0 ( mod p )

If we keep factorizing 2 r - 2 , eventaully it would become 2 1 , which would give us the below equation

=> ( a 2 1 . d - 1 ) … ( a 2 r - 2 . d + 1 ) . ( a 2 r - 1 . d + 1 ) = 0 ( mod p )

=> ( a d - 1 ) . ( a d + 1 ) … ( a 2 r - 2 . d + 1 ) . ( a 2 r - 1 . d + 1 ) = 0 ( mod p )

So this equation tells us that the number p ( which we are assuming to be prime ) has to divide atleast one of the terms on the L.H.S. But if p is not prime, then it need not divide a term on the L.H.S

Thus Miller-Rabin’s primality test check the below two conditions

• ( a d - 1 ) = 0 ( mod p ) or ( a d ) = 1 ( mod p )
• For some 0 <= t <= r - 1, ( a 2 t d + 1 ) = 0 ( mod p ) or ( a 2 t d ) = -1 ( mod p )

Example : Check if p = 221 is prime or not.
p - 1 = 220 = 2 2 . 55 . Thus we have r = 2 and d = 55.

Iteration 1 : Random base a = 174

• 174 55 mod 221 = 47. Note [ 47 != 1 ].
• When t = 0, ( 174 2 0 . 55 ) mod 221 = 47. Note [ 47 != -1 mod 221 ]
• When t = 1, ( 174 2 1 . 55 ) mod 221 = ( 174 110 ) mod 221 = 220. Note [ 220 = -1 mod 221 ] . ( Proof : -1 mod 221 = ( 220 - 221 ) mod 221 = 220 - 0 = 220 )
Thus when t = 1, it indicates that 221 could be prime or 174 is lying for 221. So try one more iteration.

Iteration 2 : Random base a = 120

• 120 55 mod 221 = 120. Note [ 120 != 1 ].
• When t = 0, ( 120 2 0 . 55 ) mod 221 = 120. Note [ 120 != -1 mod 221 ]
• When t = 1, ( 120 2 1 . 55 ) mod 221 = ( 120 110 ) mod 221 = 35. Note [ 35 != -1 mod 221 ]
Thus 120 is the witness that 221 is composite and 174 was in fact lying about 221 being prime.

Note: If we have to find 4 mod 5, we do the below
4 mod 5 = 4
4 mod 5 = ( 5 - 1 ) mod 5 = -1 mod 5

The above could be written as 4 mod 5 = -1 mod 5
So we could generalize,
If ( a ) mod p = - 1 mod p, then ( a ) mod p = ( p - 1 ) mod p.

From the above deductions we can come up with the below algorithm

Miller-Rabin’s (single test) algorithm

Check_If_Prime ( p )

1. If p equals 2 or p equals 3 then p is prime.
2. If p is even then p is definitely not prime.
3. Else find r so that p - 1 = 2 r. d for some odd d.
4. Pick a random a from the set [ 1, …, p - 1 ]
5. If ( a d ) = 1 ( mod p ) then n could be prime.
6. Else for some 0 <= t <= r - 1, if ( a 2 t d ) = -1 ( mod p ) then p could be prime.
7. Else p is composite.

Example to get an idea of the probability of getting the correct answer:
Consider p = 15, then for a = 1, 2, 3, . . . , 14
We get the values of a p − 1 - 1 % p as = [ 0, 3, 8, 0, 9, 5, 3, 3, 5, 9, 0, 8, 3, 0 ]

Thus for a ∈ { 1, 4, 11, 14 } we have a p - 1 - 1 = 0 ( mod p )
So only for 4 out of 14 numbers conform a p - 1 - 1 = 0 ( mod p ).

Conclusion :

• If p divides at least one of the terms on the L.H.S then it is probably a prime number with a probability of being prime as [ 3 / 4 ] for a single test
i.e probability of being composite [ 1 / 4 ]
.
• If p divides not a single terms on the L.H.S. it is has to be composite.

If we choose another random value for a from 1 … ( p - 1 ) the probability of Miller-Rabin’s test going wrong would go down from ( 1 / 4 ) to ( 1 / 4 ) * ( 1 / 4 ) i.e ( 1 / 16 ). Thus if we use 10 values of a the chance of going wrong would decrease from ( 1 / 4 ) to ( 1 / 4 ) 10

Implementation of Miller-Rabin’s primality test

``````import random

# p is the number to check, a is the base, d is odd, r is the power
def Check_If_Composite (p, a, d, r) :

# If the below is true, then the number is likely to be prime.
# Condition 1: If ( a ^ d ) = 1 ( mod p ) i.e find out if ( a ^ d ) % p = 1
# Condition 2: If ( a ^ ( ( 2 ^ 0 ) . d ) = - 1 ( mod p )
# i.e find out if ( a ^ d ) = -1 (mod p) which is same as finding out if ( a ^ d ) % p = p - 1
remainder = pow (a, d, p)

if (remainder == 1 or remainder == p - 1) :
return False

# Note : Remainder is already calulated above. It is (a ^ d) % p.
# Below loop would calulate if (a ^ (( 2 ^ t ).d)) == -1 (mod p) for some 1 <= t <= r-1
# And ( a ^ d ) = -1 (mod p) is same as ( a ^ d )  % p = p - 1
# Example ( a ^ d . a ^ d) % p = (a ^ 2.d) mod p
for t in range (1, r) :
if ((remainder * remainder) % p == p - 1) :
return False
return True

def Check_If_Prime (p, iterations) :

if (p == 2 or p == 3) :
return True

# p is even
if (p % 2 == 0) :
return False

# p - 1 has to be even as p is odd
d = p - 1
r = 0
while ((d % 2) == 0) :
d = int ( d / 2 )
r += 1

for i in range ( iterations + 1 ) :
# Choose a random a ∈ { 1, 2, . . . , p − 1 }.
a = random.randrange ( 2, p )
# If it is definitely composite, then it is not prime.
if ( Check_If_Composite(p, a, d, r) == True ) :
return False
return True

def main():

tests = int(input("How many numbers to check : "))

for i in range(tests):
p = int(input("Check if prime : "))
if (Check_If_Prime(p, 5) == True):
print("Yes")
else:
print("No")

if __name__ == "__main__" :
main()
``````

Output

``````How many numbers to check : 10
Check if prime : 13
Yes
Check if prime : 221
No
Check if prime : 743
Yes
Check if prime : 751
Yes
Check if prime : 1067
No
Check if prime : 1069
Yes
Check if prime : 3571
Yes
Check if prime : 3251
Yes
Check if prime : 3409
No
Check if prime : 741
No
``````
``````#include<iostream>

using namespace std;
typedef unsigned long long ULL;

// Evaluates ( a * b ) mod p. This is needed when a * b exceeds the sizeof of ULL;
ULL Modular_Mul (ULL a, ULL b, ULL p) {

ULL result = 0;
a = a % p;
while ( b > 0 ) {
if ( b & 1 ) {
result = ( result + a ) % p;
}
a = ( 2 * a ) % p;
b /= 2;
}
return result % p;
}

// Below function evaluates ( base ^ power ) mod p
ULL Modular_Expo (ULL base, ULL power, ULL p) {

ULL result = 1;
while ( power > 0 ) {
if ( power & 1 ) {
// Instead of evaluating the result as result = ( result * base ) % p
// do the multiplication under modulo
result = Modular_Mul(result, base, p);
}

// Instead of evaluating base as base = ( base * base ) % p
// do the multiplication of base under modulo
base = Modular_Mul(base, base, p);
power >>= 1; // power = power / 2
}
return result % p;
}

// p is the number to check, a is the base, d is odd, r is the power
bool Check_If_Composite (ULL p, ULL a, ULL d, ULL r) {

// If the below is true, then the number is likely to be prime.
// Condition 1: If ( a ^ d ) = 1 ( mod p ) i.e find out if ( a ^ d ) % p = 1
// Condition 2: If ( a ^ ( ( 2 ^ 0 ) . d ) = - 1 ( mod p )
// i.e find out if ( a ^ d ) = -1 (mod p) which is same as finding out if ( a ^ d ) % p = p - 1
ULL remainder = Modular_Expo(a, d, p);

if (remainder == 1 or remainder == p - 1)
return false;

// Note : Remainder is already calulated above. It is (a ^ d) % p.
// Below loop would calulate if (a ^ (( 2 ^ t ).d)) == -1 (mod p) for some 1 <= t <= r-1
// And ( a ^ d ) = -1 (mod p) is same as ( a ^ d )  % p = p - 1
// Example ( a ^ d . a ^ d) % p = (a ^ 2.d) mod p */
for (int t = 1; t < r; t++) {
remainder = Modular_Mul (remainder, remainder, p);
if (remainder == p - 1) {
return false;
}
}
return true;
}

bool Check_If_Prime (ULL p, int iterations) {

if (p == 2 or p == 3)
return true;

// p is even
if (p % 2 == 0)
return false;

// p - 1 has to be even as p is odd
ULL d = p - 1;
ULL r = 0;
while ( (d & 1) == 0 ) {
d >>= 1;
r++;
}

for (int i = 0; i <= iterations; i++) {
// Choose a random a ∈ { 1, 2, . . . , p − 1 }.
ULL a = 2 + rand() % (p - 2);
// If it is definitely composite, then it is not prime.
if ( Check_If_Composite(p, a, d, r) == true ) {
return false;
}
}
return true;
}

int main() {

int tests;
cout << "How many numbers to check: ";
cin >> tests;

ULL p;
for (int i = 0; i < tests; i++) {
cout << "Check : ";
cin >> p;
if (Check_If_Prime(p, 5))
cout << "Yes" << endl;
else
cout << "No" << endl;
}
return 0;
}
``````

Output

``````How many numbers to check: 10
Check : 13
Yes
Check : 221
No
Check : 741
No
Check : 743
Yes
Check : 751
Yes
Check : 1067
No
Check : 1069
Yes
Check : 3571
Yes
Check : 3251
Yes
Check : 3407
Yes
``````
``````import java.util.Scanner;
import java.util.Random;

class CheckPrime {
// Evaluates ( a * b ) mod p. This is needed when a * b exceeds the sizeof of int;
int Modular_Mul (int a, int b, int p) {

int result = 0;
a = a % p;
while ( b > 0 ) {
if ( b % 2 == 1 ) {
result = ( result + a ) % p;
}
a = ( 2 * a ) % p;
b /= 2;
}
return result % p;
}

// Below function evaluates ( base ^ power ) mod p
int Modular_Expo (int base, int power, int p) {

int result = 1;
while ( power > 0 ) {
if ( power % 2 == 1 ) {
// Instead of evaluating the result as result = ( result * base ) % p
// do the multiplication under modulo
result = Modular_Mul(result, base, p);
}

// Instead of evaluating base as base = ( base * base ) % p
// do the multiplication of base under modulo
base = Modular_Mul(base, base, p);
power >>= 1; // power = power / 2
}
return result % p;
}

// p is the number to check, a is the base, d is odd, r is the power
boolean Check_If_Composite (int p, int a, int d, int r) {

// If the below is true, then the number is likely to be prime.
// Condition 1: If ( a ^ d ) = 1 ( mod p ) i.e find out if ( a ^ d ) % p = 1
// Condition 2: If ( a ^ ( ( 2 ^ 0 ) . d ) = - 1 ( mod p )
// i.e find out if ( a ^ d ) = -1 (mod p) which is same as finding out if ( a ^ d ) % p = p - 1
int remainder = Modular_Expo (a, d, p);

if (remainder == 1 || remainder == p - 1)
return false;

// Note : Remainder is already calulated above. It is (a ^ d) % p.
// Below loop would calulate if (a ^ (( 2 ^ t ).d)) == -1 (mod p) for some 1 <= t <= r-1
// And ( a ^ d ) = -1 (mod p) is same as ( a ^ d )  % p = p - 1
// Example ( a ^ d . a ^ d) % p = (a ^ 2.d) mod p */
for (int t = 1; t < r; t++) {
remainder = Modular_Mul (remainder, remainder, p);
if (remainder == p - 1) {
return false;
}
}
return true;
}

boolean Check_If_Prime (int p, int iterations) {

if (p == 2 || p == 3)
return true;

// p is even
if (p % 2 == 0)
return false;

// p - 1 has to be even as p is odd
int d = p - 1;
int r = 0;
while ( (d & 1) == 0 ) {
d >>= 1;
r++;
}

Random rand = new Random();

for (int i = 0; i <= iterations; i++) {
// Choose a random a ∈ { 1, 2, . . . , p − 1 }.
int a = rand.nextInt ( p - 1 ) + 2;
// If it is definitely composite, then it is not prime.
if ( Check_If_Composite(p, a, d, r) ) {
return false;
}
}
return true;
}

public static void main (String[] args) {

Scanner sc= new Scanner(System.in);
System.out.print("How many numbers to check: ");
int tests = sc.nextInt();

CheckPrime obj = new CheckPrime();
int p;
for (int i = 0; i < tests; i++) {
System.out.print("Check if prime : ");
p = sc.nextInt();
if (obj.Check_If_Prime(p, 5))
System.out.println("Yes");
else
System.out.println("No");
}
}
}
``````

Output

``````How many numbers to check: 10
Check if prime : 13
Yes
Check if prime : 221
No
Check if prime : 741
No
Check if prime : 743
Yes
Check if prime : 751
Yes
Check if prime : 1067
No
Check if prime : 1069
Yes
Check if prime : 3571
Yes
Check if prime : 3251
Yes
Check if prime : 3407
Yes
``````