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 false 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


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



My Personal Notes arrow_drop_up

Improved By : jit_t, kartikay101