Miller–Rabin primality test: Difference between revisions

=={{header|C}}==

{{libheader|GMP}}

'''miller-rabin.h'''

<lang c>#ifndef _MILLER_RABIN_H_

#define _MILLER_RABIN_H

#include <gmp.h>

bool miller_rabin_test(mpz_t n, int j);

#endif</lang>

'''miller-rabin.c'''

{{trans|Fortran}}

For <code>decompose</code> (and header <tt>primedecompose.h</tt>), see [[Prime decomposition#C|Prime decomposition]].

<lang c>#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;

}</lang>

'''Testing'''

<lang c>#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;

}</lang>