How to resolve the algorithm Miller–Rabin primality test step by step in the C programming language

Problem Statement
Step by Step Solution
Sourcecode

Problem Statement

The Miller–Rabin primality test or Rabin–Miller primality test is a primality test: an algorithm which determines whether a given number is prime or not. The algorithm, as modified by Michael O. Rabin to avoid the generalized Riemann hypothesis, is a probabilistic algorithm. The pseudocode, from Wikipedia is:

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Miller–Rabin primality test step by step in the C programming language

The Miller-Rabin Primality Test

The Miller-Rabin primality test is a probabilistic test for primality. It is much faster than the deterministic primality tests, but it is not as reliable. The Miller-Rabin test is based on the following theorem:

If $n$ is a prime number, then for any integer $a$ in the range $1 \le a \le n-1$, the following two conditions hold:

$a^n \equiv a \pmod{n}$

If $a^d \equiv 1 \pmod{n}$ for some $d$ that divides $n-1$, then $d = n-1$.

The Algorithm

The Miller-Rabin test is a probabilistic test, meaning that it does not always return the correct answer. However, the probability of a false positive is very low. The algorithm works as follows:

Choose a random integer $a$ in the range $1 \le a \le n-1$.
Compute $a^n \bmod{n}$.
If $a^n \equiv a \pmod{n}$, then $n$ is probably prime.
If $a^n \not\equiv a \pmod{n}$, then $n$ is definitely not prime.
If $a^d \equiv 1 \pmod{n}$ for some $d$ that divides $n-1$, then $d = n-1$ and $n$ is probably prime.
If none of the above conditions hold, then $n$ is definitely not prime.

The Code

The following code implements the Miller-Rabin primality test in C:

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gmp.h>

bool miller_rabin_test(mpz_t n, int j)
{
 bool res;
 mpz_t f[MAX_DECOMPOSE];
 mpz_t s, d, a, x, r;
 mpz_t n_1, n_3;
 gmp_randstate_t rs;
 int l=0, k;

 res = false;
 gmp_randinit_default(rs);

 mpz_init(s); mpz_init(d);
 mpz_init(a); mpz_init(x); mpz_init(r);
 mpz_init(n_1); mpz_init(n_3);

 if ( mpz_cmp_si(n, 3) <= 0 ) { // let us consider 1, 2, 3 as prime
   gmp_randclear(rs);
   return true;
 }
 if ( mpz_odd_p(n) != 0 ) {
   mpz_sub_ui(n_1, n, 1);         //  n-1
   mpz_sub_ui(n_3, n, 3);         //  n-3
   l = decompose(n_1, f);
   mpz_set_ui(s, 0);
   mpz_set_ui(d, 1);
   for(k=0; k < l; k++) {
     if ( mpz_cmp_ui(f[k], 2) == 0 ) 
   mpz_add_ui(s, s, 1);
     else
   mpz_mul(d, d, f[k]);
   }                             // 2^s * d = n-1
   while(j-- > 0) {
     mpz_urandomm(a, rs, n_3);     // random from 0 to n-4
     mpz_add_ui(a, a, 2);          // random from 2 to n-2
     mpz_powm(x, a, d, n);
     if ( mpz_cmp_ui(x, 1) == 0 ) continue;
     mpz_set_ui(r, 0);
     while( mpz_cmp(r, s) < 0 ) {
   if ( mpz_cmp(x, n_1) == 0 ) break;
   mpz_powm_ui(x, x, 2, n);
   mpz_add_ui(r, r, 1);
     }
     if ( mpz_cmp(x, n_1) == 0 ) continue;
     goto flush; // woops
   }
   res = true;
 }

flush:
 for(k=0; k < l; k++) mpz_clear(f[k]);
 mpz_clear(s); mpz_clear(d);
 mpz_clear(a); mpz_clear(x); mpz_clear(r);
 mpz_clear(n_1); mpz_clear(n_3);
 gmp_randclear(rs);
 return res;
}

Usage

The Miller-Rabin primality test can be used to test the primality of any integer. The following code shows how to use the test:

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gmp.h>

int main()
{
 mpz_t n;

 mpz_init(n);
 mpz_set_ui(n, 1);
 
 while ( mpz_cmp_ui(n, TOP) < 0 ) {
   if ( miller_rabin_test(n, PREC) ) {
     gmp_printf("%Zd maybe prime\n", n);
   } /*else {
     gmp_printf("%Zd not prime\n", n);
     }*/ // remove the comment iff you're interested in
         // sure non-prime.
   mpz_add_ui(n, n, 1);
 }

 mpz_clear(n);
 return EXIT_SUCCESS;
}

The output of the program is a list of all the prime numbers less than TOP.

Source code in the c programming language

#ifndef _MILLER_RABIN_H_
#define _MILLER_RABIN_H
#include <gmp.h>
bool miller_rabin_test(mpz_t n, int j);
#endif


#include <stdbool.h>
#include <gmp.h>
#include "primedecompose.h"

#define MAX_DECOMPOSE 100

bool miller_rabin_test(mpz_t n, int j)
{
  bool res;
  mpz_t f[MAX_DECOMPOSE];
  mpz_t s, d, a, x, r;
  mpz_t n_1, n_3;
  gmp_randstate_t rs;
  int l=0, k;

  res = false;
  gmp_randinit_default(rs);

  mpz_init(s); mpz_init(d);
  mpz_init(a); mpz_init(x); mpz_init(r);
  mpz_init(n_1); mpz_init(n_3);

  if ( mpz_cmp_si(n, 3) <= 0 ) { // let us consider 1, 2, 3 as prime
    gmp_randclear(rs);
    return true;
  }
  if ( mpz_odd_p(n) != 0 ) {
    mpz_sub_ui(n_1, n, 1);         //  n-1
    mpz_sub_ui(n_3, n, 3);         //  n-3
    l = decompose(n_1, f);
    mpz_set_ui(s, 0);
    mpz_set_ui(d, 1);
    for(k=0; k < l; k++) {
      if ( mpz_cmp_ui(f[k], 2) == 0 ) 
	mpz_add_ui(s, s, 1);
      else
	mpz_mul(d, d, f[k]);
    }                             // 2^s * d = n-1
    while(j-- > 0) {
      mpz_urandomm(a, rs, n_3);     // random from 0 to n-4
      mpz_add_ui(a, a, 2);          // random from 2 to n-2
      mpz_powm(x, a, d, n);
      if ( mpz_cmp_ui(x, 1) == 0 ) continue;
      mpz_set_ui(r, 0);
      while( mpz_cmp(r, s) < 0 ) {
	if ( mpz_cmp(x, n_1) == 0 ) break;
	mpz_powm_ui(x, x, 2, n);
	mpz_add_ui(r, r, 1);
      }
      if ( mpz_cmp(x, n_1) == 0 ) continue;
      goto flush; // woops
    }
    res = true;
  }

flush:
  for(k=0; k < l; k++) mpz_clear(f[k]);
  mpz_clear(s); mpz_clear(d);
  mpz_clear(a); mpz_clear(x); mpz_clear(r);
  mpz_clear(n_1); mpz_clear(n_3);
  gmp_randclear(rs);
  return res;
}


#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <gmp.h>
#include "miller-rabin.h"

#define PREC 10
#define TOP  4000

int main()
{
  mpz_t num;

  mpz_init(num);
  mpz_set_ui(num, 1);
  
  while ( mpz_cmp_ui(num, TOP) < 0 ) {
    if ( miller_rabin_test(num, PREC) ) {
      gmp_printf("%Zd maybe prime\n", num);
    } /*else {
      gmp_printf("%Zd not prime\n", num);
      }*/ // remove the comment iff you're interested in
          // sure non-prime.
    mpz_add_ui(num, num, 1);
  }

  mpz_clear(num);
  return EXIT_SUCCESS;
}


// calcul a^n%mod
size_t power(size_t a, size_t n, size_t mod)
{
    size_t power = a;
    size_t result = 1;

    while (n)
    {
        if (n & 1)
            result = (result * power) % mod;
        power = (power * power) % mod;
        n >>= 1;
    }
    return result;
}

// n−1 = 2^s * d with d odd by factoring powers of 2 from n−1
bool witness(size_t n, size_t s, size_t d, size_t a)
{
    size_t x = power(a, d, n);
    size_t y;

    while (s) {
        y = (x * x) % n;
        if (y == 1 && x != 1 && x != n-1)
            return false;
        x = y;
        --s;
    }
    if (y != 1)
        return false;
    return true;
}

/*
 * if n < 1,373,653, it is enough to test a = 2 and 3;
 * if n < 9,080,191, it is enough to test a = 31 and 73;
 * if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
 * if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;
 * if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
 * if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
 * if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.
 */

bool is_prime_mr(size_t n)
{
    if (((!(n & 1)) && n != 2 ) || (n < 2) || (n % 3 == 0 && n != 3))
        return false;
    if (n <= 3)
        return true;

    size_t d = n / 2;
    size_t s = 1;
    while (!(d & 1)) {
        d /= 2;
        ++s;
    }

    if (n < 1373653)
        return witness(n, s, d, 2) && witness(n, s, d, 3);
    if (n < 9080191)
        return witness(n, s, d, 31) && witness(n, s, d, 73);
    if (n < 4759123141)
        return witness(n, s, d, 2) && witness(n, s, d, 7) && witness(n, s, d, 61);
    if (n < 1122004669633)
        return witness(n, s, d, 2) && witness(n, s, d, 13) && witness(n, s, d, 23) && witness(n, s, d, 1662803);
    if (n < 2152302898747)
        return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11);
    if (n < 3474749660383)
        return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13);
    return witness(n, s, d, 2) && witness(n, s, d, 3) && witness(n, s, d, 5) && witness(n, s, d, 7) && witness(n, s, d, 11) && witness(n, s, d, 13) && witness(n, s, d, 17);
}


typedef unsigned long long int ulong;

ulong mul_mod(ulong a, ulong b, const ulong mod) {
	ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
	for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
	return res % mod;
}

ulong pow_mod(ulong n, ulong exp, const ulong mod) {
	ulong res = 1; // return (n ^ exp) % mod
	for (n %= mod; exp; exp & 1 ? res = mul_mod(res, n, mod) : 0, n = mul_mod(n, n, mod), exp >>= 1);
	return res;
}

int is_prime(ulong N) {
	// Perform a Miller-Rabin test, it should be a deterministic version.
	const ulong n_primes = 9, primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
	for (ulong i = 0; i < n_primes; ++i)
		if (N % primes[i] == 0) return N == primes[i];
	if (N < primes[n_primes - 1]) return 0;
	int res = 1, s = 0;
	ulong t;
	for (t = N - 1; ~t & 1; t >>= 1, ++s);
	for (ulong i = 0; i < n_primes && res; ++i) {
		ulong B = pow_mod(primes[i], t, N);
		if (B != 1) {
			for (int b = s; b-- && (res = B + 1 != N);)
				B = mul_mod(B, B, N);
			res = !res;
		}
	}
	return res;
}

int main(void){
	return is_prime(8193145868754512737);
}

You may also check:How to resolve the algorithm Increment a numerical string step by step in the Python programming language
You may also check:How to resolve the algorithm The Twelve Days of Christmas step by step in the Jsish programming language
You may also check:How to resolve the algorithm HTTPS/Authenticated step by step in the Raku programming language
You may also check:How to resolve the algorithm Sort stability step by step in the Pascal programming language
You may also check:How to resolve the algorithm Generate lower case ASCII alphabet step by step in the AutoHotkey programming language

How to resolve the algorithm Miller–Rabin primality test step by step in the C programming language

Table of Contents

Problem Statement

Step by Step solution about How to resolve the algorithm Miller–Rabin primality test step by step in the C programming language

Source code in the c programming language