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.

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

 // C++ program Miller-Rabin primality test #include 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*2r = 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: ";     for (int n = 1; n < 100; n++)        if (isPrime(n, k))           cout << n << " ";        return 0; }

/div>

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*2r     // = 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.*/

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*2r = 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

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*2r     // = 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

PHP

 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*2r = 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: ";     for (\$n = 1; \$n < 100; \$n++)     if (isPrime(\$n, \$k))         echo \$n , " ";    // This code is contributed by ajit ?>

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
OR
for some i, where 0 <= i <= r-1.

Next Article :
Primality Test | Set 4 (Solovay-Strassen)