Primality Test | Set 3 (Miller–Rabin)

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, Miller-Rabin method is discussed. This method is a probabilistic method (Like Fermat), but it 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 above algorithm.

C++

filter_none

edit
close

play_arrow

link
brightness_4
code

// 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 false if n is
// probably prime.
// d is an odd number such that  d*2<sup>r</sup> = 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 nber 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;
}

chevron_right


Java

filter_none

edit
close

play_arrow

link
brightness_4
code

// 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 nber 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.*/

chevron_right


Python3

filter_none

edit
close

play_arrow

link
brightness_4
code

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

chevron_right


C#

filter_none

edit
close

play_arrow

link
brightness_4
code

// 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 nber 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

chevron_right


PHP

filter_none

edit
close

play_arrow

link
brightness_4
code

<?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 nber 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
?>

chevron_right



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 

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 above two points, for every randomly picked number in range [2, n-2], 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)

References:
https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

This article is contributed Ruchir Garg. Please write comments if you find anything incorrect, or you want to share more information about the topic discussed above

GeeksforGeeks has prepared a complete interview preparation course with premium videos, theory, practice problems, TA support and many more features. Please refer Placement 100 for details




My Personal Notes arrow_drop_up