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.

There are some composite numbers that also pass the Fermat’s test a p - 1 - 1 = 0 ( mod p ) for all values of a i.e a ∈ { 1, 2, . . . , p − 1 } .
Such composite numbers are called Carmichael numbers.
The first seven Carmichael numbers are 561, 1105, 1729, 2465, 2821, 6601, and 8911; and there are infinitely many such Carmichael numbers.

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

Copyright (c) 2019-2023, Algotree.org.
All rights reserved.