Miller-Rabin primality test (Groovy)

From LiteratePrograms

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 / 4^k 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.

We will implement the test for arbitrary-precision integers, with Java's BigInteger type, which is also used more easily in Groovy.

When performing the Miller-Rabin pass, only a, n, d, and powers of a are large numbers; s and i are logarithmically smaller. We precompute the large integers 1 and n−1 for use in comparisons and squaring:

<<Miller-Rabin pass>>=
def miller_rabin_pass(BigInteger a, BigInteger n) {
    BigInteger d, s, a_to_power
    compute s and d
    a_to_power = a.modPow(d, n)
    if (a_to_power == 1) {
        return true
    }
    for (i in 0..<s) {
        if (a_to_power == n - 1) {
            return true
        }
        a_to_power = (a_to_power * a_to_power) % n
    }
    return (a_to_power == n - 1)
}

We've replaced the squaring using modular_exponent with a modular multiply of a_to_power by itself. There does not exist the right shift operation on the Java's BigInteger type. So we have used d.divide(2G), where 2G is Groovy's literal for the BigInteger constant 2.

<<compute s and d>>=
d = n - 1
s = 0
while (d % 2 == 0) {
    d = d.divide(2G)
    s += 1
}

Since our numbers are now arbitrarily large, we can no longer take advantage of Jaeschke's result to produce a deterministic implementation, but we can produce a highly accurate probabilistic implementation by simply running the test enough times to make 1/4^k very small, say 20 times.

<<Miller-Rabin>>=
def miller_rabin(n) {
    def k = 20
    for (i in 0..<k) {
        BigInteger a = 0
        Random rand = null
        while (a == 0) {
            rand = new Random(System.currentTimeMillis())
             a = new BigDecimal(rand.nextDouble()).multiply(n).toBigInteger()
        }
        if (!miller_rabin_pass(a, n)) {
            return false
        }
    }
    return true
}

For convenience we ignore the special cases where n ≤ 2.

Finally, some test code:

<<MillerRabin.groovy>>=
Miller-Rabin pass
Miller-Rabin
BigInteger n, nbits, p
if (args[0] == "test") {
    n = new BigInteger(args[1])
    println (miller_rabin(n) ? "PRIME" : "COMPOSITE")
}
else if (args[0] == "genprime") {
    nbits = new BigInteger(args[1])
    Random rand = null
    while (true) {
        rand = new Random(System.currentTimeMillis())
        p = new BigInteger(nbits.intValue(), rand)
        if (p % 2 == 0)
            continue
	else if (p % 3 == 0)
            continue
	else if (p % 5 == 0)
            continue
	else if (p % 7 == 0)
            continue
        if (miller_rabin(p)) {
            println p
            break
        }
    }
}

We've augmented main() with the ability to generate a prime number of a specified number of bits. It does this by randomly selecting numbers of that size 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 (p % 2 == 0)
            continue
	else if (p % 3 == 0)
            continue
	else if (p % 5 == 0)
            continue
	else if (p % 7 == 0)
            continue

Here's some sample output:

$ groovy MillerRabin.groovy test 516119616549881
PRIME
$ groovy MillerRabin.groovy test 516119616549887
COMPOSITE
$ groovy MillerRabin.groovy genprime 128
338012683331153239276905527729082198073
$ groovy MillerRabin.groovy genprime 512
79341348271742693405644860279850296923552898454148485677160734904706128744200302\
42901825660196095090061259225626245698999291562967073071438566777530365841
$ groovy MillerRabin.groovy genprime 1500
59579073968246513417478073620871442204747950333852855582769347422962540842369159\
46073405832930527012381929952727658841800423072204289685129481530589571531880578\
91100103693890987958207702581355160770537505615147869598330091266089878569072355\
72228424009915737779322260330614377520077501154638734354238083483760688090842520\
12237118646930534600916857916698772040032336635790847660021841905365318564285442\
998224235824474092228111652583101050232658251560933

Download code