Miller-Rabin primality test (C, GMP)

From LiteratePrograms

Jump to: navigation, search
Other implementations: C | C, GMP | Clojure | Groovy | Java | Python | Ruby | Scala

The Miller-Rabin primality test is a simple probabilistic algorithm for determining whether a number is prime or composite that is easy to implement. It proves compositeness of a number using the following formulas:

Suppose 0 < a < n is coprime to n (this is easy to test using the GCD). Write the number n−1 as , where d is odd. Then, provided that all of the following formulas hold, n is composite:

for all

If a is chosen uniformly at random and n is prime, these formulas hold with probability 1/4. Thus, repeating the test for k random choices of a gives a probability of 1 − 1 / 4k that the number is prime. Moreover, Gerhard Jaeschke showed that any 32-bit number can be deterministically tested for primality by trying only a=2, 7, and 61.


This C implementation uses the GNU Multiple Precision (GMP) arithmetic library. Using GMP, the test can be implemented in an efficient and relatively straightforward manner based on the definition. Note that GMP already supplies a good primality testing function, mpz_probab_prime_p, which is based on Miller-Rabin (as described in the manual); the supplied function should be preferred, and the following code is only for illustration.

We write a function that tests the formulas above for a given a and n. It works by starting with ad and squaring it s−1 times, comparing each to −1 (n−1 actually; we have to add n because we're working with nonnegative numbers). To compute the power ad quickly we take advantage of GMP's modular exponentiation function mpz_powm:

<<Miller-Rabin pass>>=
int miller_rabin_pass(mpz_t a, mpz_t n) {
    int i, s, result;
    mpz_t a_to_power, d, n_minus_one;
    mpz_init(n_minus_one);
    mpz_sub_ui(n_minus_one, n, 1);
    compute s and d
    mpz_init(a_to_power);
    mpz_powm(a_to_power, a, d, n);
    if (mpz_cmp_ui(a_to_power, 1) == 0)  {
        result=PROBABLE_PRIME; goto exit;
    }
    for(i=0; i < s-1; i++) {
        if (mpz_cmp(a_to_power, n_minus_one) == 0) {
            result=PROBABLE_PRIME; goto exit;
        }
        mpz_powm_ui (a_to_power, a_to_power, 2, n);
    }
    if (mpz_cmp(a_to_power, n_minus_one) == 0) {
        result=PROBABLE_PRIME; goto exit;
    }
    result = COMPOSITE;
exit:
    mpz_clear(a_to_power);
    mpz_clear(d);
    mpz_clear(n_minus_one);
    return result;
}

Above, only a, n, d, and powers of a are large numbers; s and i are logarithmically smaller and so can be represented with machine words. We precompute the large integer n−1 for use in comparisons. The return values are the same as those used by mpz_probab_prime_p:

<<Miller-Rabin return values>>=
#define COMPOSITE        0
#define PROBABLE_PRIME   1

To compute s and d, we just divide n−1 by 2 until it becomes odd. We take advantage of mpz_fdiv_q_2exp for efficient division by 2:

<<compute s and d>>=
s = 0;
mpz_init_set(d, n_minus_one);
while (mpz_even_p(d)) {
    mpz_fdiv_q_2exp(d, d, 1);
    s++;
}

A single pass has, in theory, up to a 1/4 chance of being incorrect. We write a wrapper function that calls it a number of times with different random a to shrink this probability to something very tiny. We use GMP's mpz_urandomm to produce the random a, rejecting any zeros. Since we want to use a difference random sequence for each call, we require the random state to be passed in.

<<Miller-Rabin>>=
int miller_rabin(mpz_t n, gmp_randstate_t rand_state) {
    mpz_t a;
    int repeat;
    mpz_init(a);
    for(repeat=0; repeat<20; repeat++) {
        do {
            mpz_urandomm(a, rand_state, n);
        } while (mpz_sgn(a) == 0);
        if (miller_rabin_pass(a, n) == COMPOSITE) {
            return COMPOSITE;
        }
    }
    return PROBABLE_PRIME;
}

Since small n ≤ 2 would normally be tested with trial division anyway, we ignore these special cases.

Finally, some test code:

<<miller-rabin.c>>=
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include <time.h>
Miller-Rabin return values
Miller-Rabin pass
Miller-Rabin
int main(int argc, char* argv[]) {
    mpz_t n, max, two, p;
    gmp_randstate_t rand_state;
    gmp_randinit_default(rand_state);
    gmp_randseed_ui (rand_state, time(NULL));
    if (strcmp(argv[1], "test") == 0) {
        mpz_init_set_str(n, argv[2], 10);
        puts(miller_rabin(n, rand_state) == PROBABLE_PRIME ? "PRIME" : "COMPOSITE");
    } else if (strcmp(argv[1], "genprime") == 0) {
        mpz_init(max);
        mpz_init_set_ui(two, 2);
        mpz_mul_2exp(max, two, atoi(argv[2])+1);
        mpz_init(p);
        do {
            mpz_urandomm(p, rand_state, max);
            test for small factors
        } while (miller_rabin(p, rand_state) == COMPOSITE);
        mpz_out_str(stdout, 10, p);
        puts("");
    }
    return 0;
}

This test main() performs two main tasks: it can test a given number for primality, or it can generate a prime of a specified number of bits by choosing random numbers until it finds a prime one. We could use this to rapidly generate large primes for use in cryptography. Since most random values have small prime factors, we first test some of these to avoid an expensive test:

<<test for small factors>>=
if (mpz_even_p(p)) continue;
if (mpz_fdiv_ui(p, 3) == 0) continue;
if (mpz_fdiv_ui(p, 5) == 0) continue;
if (mpz_fdiv_ui(p, 7) == 0) continue;

Here's some sample output using the default random seed:

$ ./miller-rabin test 516119616549881
PRIME
$ ./miller-rabin test 516119616549887
COMPOSITE
$ ./miller-rabin genprime 128
750939670750127524461407364479090229319
$ ./miller-rabin genprime 512
3197120138436173701226583134640728580338949002968072415067612719491901934327385\
213861339206923061932145554021154341263829908046193368844987237491712679531
$ ./miller-rabin genprime 1500
7774170786735403788138484057436891511278248835103163926563116623429983126719713\
2508562677357193939022326689674599429496751069048535076610006998558461169350187\
0213636530383574519176090749776752244752659815375020514536232328587844581509305\
6530497921265833050905987199183587066829302220365425705920279080352809428940203\
0929245116597109258323265772559335793346530534939803881784890423036766501992454\
57216718122569053613650716247455521949998232588949376723

Independent confirmation shows these results are correct. On a Pentium 4 2.8 GhZ machine, the last run took between 14 and 90 seconds. Runtime varies widely from run to run depending on how quickly we stumble across a prime; by the prime number theorem, on average about 450 random numbers have to be tested, of which about 70% are rejected immediately by the trial division tests. Using GMP's much faster built-in mpz_probab_prime_p, the same run takes between 1/2 and 6 seconds; although the algorithm is the same, it is finely tuned for performance.

Both implementations are much faster than the one in Miller-Rabin primality test (C) based on very simplistic arbitrary-precision arithmetic. Even the slower of the two takes only 0.2 to 3 seconds to generate a 512-bit prime, whereas the implementation not based on GMP required in excess of 20 seconds.

Download code
Views