Open In App

Primality Test | Set 3 (Miller–Rabin)

Improve
Improve
Improve
Like Article
Like
Save Article
Save
Share
Report issue
Report

Given a number n, check if it is prime or not. We have introduced and discussed School and Fermat methods for primality testing.
Primality Test | Set 1 (Introduction and School Method) 
Primality Test | Set 2 (Fermat Method)
In this post, the Miller-Rabin method is discussed. This method is a probabilistic method ( like Fermat), but it is generally preferred over Fermat’s method.

Algorithm:

// It returns false if n is composite and returns true if n
// is probably prime.  k is an input parameter that determines
// accuracy level. Higher value of k indicates more accuracy.
bool isPrime(int n, int k)
1) Handle base cases for n < 3
2) If n is even, return false.
3) Find an odd number d such that n-1 can be written as d*2r. 
   Note that since n is odd, (n-1) must be even and r must be 
   greater than 0.
4) Do following k times
     if (millerTest(n, d) == false)
          return false
5) Return true.

// This function is called for all k trials. It returns 
// false if n is composite and returns true if n is probably
// prime.  
// d is an odd number such that d*2r = n-1 for some r>=1
bool millerTest(int n, int d)
1) Pick a random number 'a' in range [2, n-2]
2) Compute: x = pow(a, d) % n
3) If x == 1 or x == n-1, return true.

// Below loop mainly runs 'r-1' times.
4) Do following while d doesn't become n-1.
     a) x = (x*x) % n.
     b) If (x == 1) return false.
     c) If (x == n-1) return true. 

Example: 

Input: n = 13,  k = 2.

1) Compute d and r such that d*2r = n-1, 
     d = 3, r = 2. 
2) Call millerTest k times.

1st Iteration:
1) Pick a random number 'a' in range [2, n-2]
      Suppose a = 4

2) Compute: x = pow(a, d) % n
     x = 43 % 13 = 12

3) Since x = (n-1), return true.

IInd Iteration:
1) Pick a random number 'a' in range [2, n-2]
      Suppose a = 5

2) Compute: x = pow(a, d) % n
     x = 53 % 13 = 8

3) x neither 1 nor 12.

4) Do following (r-1) = 1 times
   a) x = (x * x) % 13 = (8 * 8) % 13 = 12
   b) Since x = (n-1), return true.

Since both iterations return true, we return true. 

Implementation: 
Below is the implementation of the above algorithm. 

C++




// C++ program Miller-Rabin primality test
#include <bits/stdc++.h>
using namespace std;
 
// Utility function to do modular exponentiation.
// It returns (x^y) % p
int power(int x, unsigned int y, int p)
{
    int res = 1;      // Initialize result
    x = x % p;  // Update x if it is more than or
                // equal to p
    while (y > 0)
    {
        // If y is odd, multiply x with result
        if (y & 1)
            res = (res*x) % p;
 
        // y must be even now
        y = y>>1; // y = y/2
        x = (x*x) % p;
    }
    return res;
}
 
// This function is called for all k trials. It returns
// false if n is composite and returns true if n is
// probably prime.
// d is an odd number such that  d*2 = n-1
// for some r >= 1
bool miillerTest(int d, int n)
{
    // Pick a random number in [2..n-2]
    // Corner cases make sure that n > 4
    int a = 2 + rand() % (n - 4);
 
    // Compute a^d % n
    int x = power(a, d, n);
 
    if (x == 1  || x == n-1)
       return true;
 
    // Keep squaring x while one of the following doesn't
    // happen
    // (i)   d does not reach n-1
    // (ii)  (x^2) % n is not 1
    // (iii) (x^2) % n is not n-1
    while (d != n-1)
    {
        x = (x * x) % n;
        d *= 2;
 
        if (x == 1)      return false;
        if (x == n-1)    return true;
    }
 
    // Return composite
    return false;
}
 
// It returns false if n is composite and returns true if n
// is probably prime.  k is an input parameter that determines
// accuracy level. Higher value of k indicates more accuracy.
bool isPrime(int n, int k)
{
    // Corner cases
    if (n <= 1 || n == 4)  return false;
    if (n <= 3) return true;
 
    // Find r such that n = 2^d * r + 1 for some r >= 1
    int d = n - 1;
    while (d % 2 == 0)
        d /= 2;
 
    // Iterate given number of 'k' times
    for (int i = 0; i < k; i++)
         if (!miillerTest(d, n))
              return false;
 
    return true;
}
 
// Driver program
int main()
{
    int k = 4;  // Number of iterations
 
    cout << "All primes smaller than 100: \n";
    for (int n = 1; n < 100; n++)
       if (isPrime(n, k))
          cout << n << " ";
 
    return 0;
}


Java




// Java program Miller-Rabin primality test
import java.io.*;
import java.math.*;
 
class GFG {
 
    // Utility function to do modular
    // exponentiation. It returns (x^y) % p
    static int power(int x, int y, int p) {
         
        int res = 1; // Initialize result
         
        //Update x if it is more than or
        // equal to p
        x = x % p;
 
        while (y > 0) {
             
            // If y is odd, multiply x with result
            if ((y & 1) == 1)
                res = (res * x) % p;
         
            // y must be even now
            y = y >> 1; // y = y/2
            x = (x * x) % p;
        }
         
        return res;
    }
     
    // This function is called for all k trials.
    // It returns false if n is composite and
    // returns false if n is probably prime.
    // d is an odd number such that d*2<sup>r</sup>
    // = n-1 for some r >= 1
    static boolean miillerTest(int d, int n) {
         
        // Pick a random number in [2..n-2]
        // Corner cases make sure that n > 4
        int a = 2 + (int)(Math.random() % (n - 4));
     
        // Compute a^d % n
        int x = power(a, d, n);
     
        if (x == 1 || x == n - 1)
            return true;
     
        // Keep squaring x while one of the
        // following doesn't happen
        // (i) d does not reach n-1
        // (ii) (x^2) % n is not 1
        // (iii) (x^2) % n is not n-1
        while (d != n - 1) {
            x = (x * x) % n;
            d *= 2;
         
            if (x == 1)
                return false;
            if (x == n - 1)
                return true;
        }
     
        // Return composite
        return false;
    }
     
    // It returns false if n is composite
    // and returns true if n is probably
    // prime. k is an input parameter that
    // determines accuracy level. Higher
    // value of k indicates more accuracy.
    static boolean isPrime(int n, int k) {
         
        // Corner cases
        if (n <= 1 || n == 4)
            return false;
        if (n <= 3)
            return true;
     
        // Find r such that n = 2^d * r + 1
        // for some r >= 1
        int d = n - 1;
         
        while (d % 2 == 0)
            d /= 2;
     
        // Iterate given number of 'k' times
        for (int i = 0; i < k; i++)
            if (!miillerTest(d, n))
                return false;
     
        return true;
    }
     
    // Driver program
    public static void main(String args[]) {
         
        int k = 4; // Number of iterations
     
        System.out.println("All primes smaller "
                                + "than 100: ");
                                 
        for (int n = 1; n < 100; n++)
            if (isPrime(n, k))
                System.out.print(n + " ");
    }
}
 
/* This code is contributed by Nikita Tiwari.*/


Python3




# Python3 program Miller-Rabin primality test
import random
 
# Utility function to do
# modular exponentiation.
# It returns (x^y) % p
def power(x, y, p):
     
    # Initialize result
    res = 1;
     
    # Update x if it is more than or
    # equal to p
    x = x % p;
    while (y > 0):
         
        # If y is odd, multiply
        # x with result
        if (y & 1):
            res = (res * x) % p;
 
        # y must be even now
        y = y>>1; # y = y/2
        x = (x * x) % p;
     
    return res;
 
# This function is called
# for all k trials. It returns
# false if n is composite and
# returns false if n is
# probably prime. d is an odd
# number such that d*2<sup>r</sup> = n-1
# for some r >= 1
def miillerTest(d, n):
     
    # Pick a random number in [2..n-2]
    # Corner cases make sure that n > 4
    a = 2 + random.randint(1, n - 4);
 
    # Compute a^d % n
    x = power(a, d, n);
 
    if (x == 1 or x == n - 1):
        return True;
 
    # Keep squaring x while one
    # of the following doesn't
    # happen
    # (i) d does not reach n-1
    # (ii) (x^2) % n is not 1
    # (iii) (x^2) % n is not n-1
    while (d != n - 1):
        x = (x * x) % n;
        d *= 2;
 
        if (x == 1):
            return False;
        if (x == n - 1):
            return True;
 
    # Return composite
    return False;
 
# It returns false if n is
# composite and returns true if n
# is probably prime. k is an
# input parameter that determines
# accuracy level. Higher value of
# k indicates more accuracy.
def isPrime( n, k):
     
    # Corner cases
    if (n <= 1 or n == 4):
        return False;
    if (n <= 3):
        return True;
 
    # Find r such that n =
    # 2^d * r + 1 for some r >= 1
    d = n - 1;
    while (d % 2 == 0):
        d //= 2;
 
    # Iterate given number of 'k' times
    for i in range(k):
        if (miillerTest(d, n) == False):
            return False;
 
    return True;
 
# Driver Code
# Number of iterations
k = 4;
 
print("All primes smaller than 100: ");
for n in range(1,100):
    if (isPrime(n, k)):
        print(n , end=" ");
 
# This code is contributed by mits


C#




// C# program Miller-Rabin primality test
using System;
 
class GFG
{
 
    // Utility function to do modular
    // exponentiation. It returns (x^y) % p
    static int power(int x, int y, int p)
    {
         
        int res = 1; // Initialize result
         
        // Update x if it is more than
        // or equal to p
        x = x % p;
 
        while (y > 0)
        {
             
            // If y is odd, multiply x with result
            if ((y & 1) == 1)
                res = (res * x) % p;
         
            // y must be even now
            y = y >> 1; // y = y/2
            x = (x * x) % p;
        }
         
        return res;
    }
     
    // This function is called for all k trials.
    // It returns false if n is composite and
    // returns false if n is probably prime.
    // d is an odd number such that d*2<sup>r</sup>
    // = n-1 for some r >= 1
    static bool miillerTest(int d, int n)
    {
         
        // Pick a random number in [2..n-2]
        // Corner cases make sure that n > 4
        Random r = new Random();
        int a = 2 + (int)(r.Next() % (n - 4));
     
        // Compute a^d % n
        int x = power(a, d, n);
     
        if (x == 1 || x == n - 1)
            return true;
     
        // Keep squaring x while one of the
        // following doesn't happen
        // (i) d does not reach n-1
        // (ii) (x^2) % n is not 1
        // (iii) (x^2) % n is not n-1
        while (d != n - 1)
        {
            x = (x * x) % n;
            d *= 2;
         
            if (x == 1)
                return false;
            if (x == n - 1)
                return true;
        }
     
        // Return composite
        return false;
    }
     
    // It returns false if n is composite
    // and returns true if n is probably
    // prime. k is an input parameter that
    // determines accuracy level. Higher
    // value of k indicates more accuracy.
    static bool isPrime(int n, int k)
    {
         
        // Corner cases
        if (n <= 1 || n == 4)
            return false;
        if (n <= 3)
            return true;
     
        // Find r such that n = 2^d * r + 1
        // for some r >= 1
        int d = n - 1;
         
        while (d % 2 == 0)
            d /= 2;
     
        // Iterate given number of 'k' times
        for (int i = 0; i < k; i++)
            if (miillerTest(d, n) == false)
                return false;
     
        return true;
    }
     
    // Driver Code
    static void Main()
    {
        int k = 4; // Number of iterations
     
        Console.WriteLine("All primes smaller " +
                                   "than 100: ");
                                 
        for (int n = 1; n < 100; n++)
            if (isPrime(n, k))
                Console.Write(n + " ");
    }
}
 
// This code is contributed by mits


PHP




<?php
// PHP program Miller-Rabin primality test
 
// Utility function to do
// modular exponentiation.
// It returns (x^y) % p
function power($x, $y, $p)
{
     
    // Initialize result
    $res = 1;
     
    // Update x if it is more than or
    // equal to p
    $x = $x % $p;
    while ($y > 0)
    {
         
        // If y is odd, multiply
        // x with result
        if ($y & 1)
            $res = ($res*$x) % $p;
 
        // y must be even now
        $y = $y>>1; // $y = $y/2
        $x = ($x*$x) % $p;
    }
    return $res;
}
 
// This function is called
// for all k trials. It returns
// false if n is composite and
// returns false if n is
// probably prime. d is an odd
// number such that d*2<sup>r</sup> = n-1
// for some r >= 1
function miillerTest($d, $n)
{
     
    // Pick a random number in [2..n-2]
    // Corner cases make sure that n > 4
    $a = 2 + rand() % ($n - 4);
 
    // Compute a^d % n
    $x = power($a, $d, $n);
 
    if ($x == 1 || $x == $n-1)
    return true;
 
    // Keep squaring x while one
    // of the following doesn't
    // happen
    // (i) d does not reach n-1
    // (ii) (x^2) % n is not 1
    // (iii) (x^2) % n is not n-1
    while ($d != $n-1)
    {
        $x = ($x * $x) % $n;
        $d *= 2;
 
        if ($x == 1)     return false;
        if ($x == $n-1) return true;
    }
 
    // Return composite
    return false;
}
 
// It returns false if n is
// composite and returns true if n
// is probably prime. k is an
// input parameter that determines
// accuracy level. Higher value of
// k indicates more accuracy.
function isPrime( $n, $k)
{
     
    // Corner cases
    if ($n <= 1 || $n == 4) return false;
    if ($n <= 3) return true;
 
    // Find r such that n =
    // 2^d * r + 1 for some r >= 1
    $d = $n - 1;
    while ($d % 2 == 0)
        $d /= 2;
 
    // Iterate given number of 'k' times
    for ($i = 0; $i < $k; $i++)
        if (!miillerTest($d, $n))
            return false;
 
    return true;
}
 
    // Driver Code
    // Number of iterations
    $k = 4;
 
    echo "All primes smaller than 100: \n";
    for ($n = 1; $n < 100; $n++)
    if (isPrime($n, $k))
        echo $n , " ";
 
// This code is contributed by ajit
?>


Javascript




<script>
// Javascript program Miller-Rabin primality test
 
// Utility function to do
// modular exponentiation.
// It returns (x^y) % p
function power(x, y, p)
{
     
    // Initialize result
    let res = 1;
     
    // Update x if it is more than or
    // equal to p
    x = x % p;
    while (y > 0)
    {
         
        // If y is odd, multiply
        // x with result
        if (y & 1)
            res = (res*x) % p;
 
        // y must be even now
        y = y>>1; // y = y/2
        x = (x*x) % p;
    }
    return res;
}
 
// This function is called
// for all k trials. It returns
// false if n is composite and
// returns false if n is
// probably prime. d is an odd
// number such that d*2<sup>r</sup> = n-1
// for some r >= 1
function miillerTest(d, n)
{
     
    // Pick a random number in [2..n-2]
    // Corner cases make sure that n > 4
    let a = 2 + Math.floor(Math.random() * (n-2)) % (n - 4);
 
    // Compute a^d % n
    let x = power(a, d, n);
 
    if (x == 1 || x == n-1)
        return true;
 
    // Keep squaring x while one
    // of the following doesn't
    // happen
    // (i) d does not reach n-1
    // (ii) (x^2) % n is not 1
    // (iii) (x^2) % n is not n-1
    while (d != n-1)
    {
        x = (x * x) % n;
        d *= 2;
 
        if (x == 1)    
            return false;
        if (x == n-1)
              return true;
    }
 
    // Return composite
    return false;
}
 
// It returns false if n is
// composite and returns true if n
// is probably prime. k is an
// input parameter that determines
// accuracy level. Higher value of
// k indicates more accuracy.
function isPrime( n, k)
{
     
    // Corner cases
    if (n <= 1 || n == 4) return false;
    if (n <= 3) return true;
 
    // Find r such that n =
    // 2^d * r + 1 for some r >= 1
    let d = n - 1;
    while (d % 2 == 0)
        d /= 2;
 
    // Iterate given number of 'k' times
    for (let i = 0; i < k; i++)
        if (!miillerTest(d, n))
            return false;
 
    return true;
}
 
// Driver Code
// Number of iterations
let k = 4;
 
document.write("All primes smaller than 100: <br>");
for (let n = 1; n < 100; n++)
    if (isPrime(n, k))
        document.write(n , " ");
 
// This code is contributed by gfgking
</script>


Output:  

All primes smaller than 100: 
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 
61 67 71 73 79 83 89 97 

Time Complexity: O(k*logn)

Auxiliary Space: O(1)

How does this work? 
Below are some important facts behind the algorithm:  

  1. Fermat’s theorem states that, If n is a prime number, then for every a, 1 <= a < n, an-1 % n = 1
  2. Base cases make sure that n must be odd. Since n is odd, n-1 must be even. And an even number can be written as d * 2s where d is an odd number and s > 0.
  3. From the above two points, for every randomly picked number in the range [2, n-2], the value of ad*2r % n must be 1.
  4. As per Euclid’s Lemma, if x2 % n = 1 or (x2 – 1) % n = 0 or (x-1)(x+1)% n = 0. Then, for n to be prime, either n divides (x-1) or n divides (x+1). Which means either x % n = 1 or x % n = -1.
  5. From points 2 and 3, we can conclude
    For n to be prime, either
    ad % n = 1 
         OR 
    ad*2i % n = -1 
    for some i, where 0 <= i <= r-1.

Next Article : 
Primality Test | Set 4 (Solovay-Strassen)
This article is contributed Ruchir Garg.
 



Last Updated : 14 Nov, 2022
Like Article
Save Article
Previous
Next
Share your thoughts in the comments
Similar Reads