Miller-Rabin primality test (Ruby)

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 Ruby's int type, which is automatically used when integers get big enough, so we don't have to do anything special.

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(a, n)
    compute s and d
    a_to_power = modPow(a, d, n)
    if a_to_power == 1
        return true
    end
    for i in 0...s do
        if a_to_power == n - 1
            return true
        end
        a_to_power = (a_to_power * a_to_power) % n
    end
    return (a_to_power == n - 1)
end

We've replaced the squaring using modular_exponent with a modular multiply of a_to_power by itself. Finally, we use bit shifts to compute d rapidly, and test just its least significant component to determine if it is odd:

<<compute s and d>>=
d = n - 1
s = 0
while d % 2 == 0 do
    d >>= 1
    s += 1
end

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)
    k = 20
    for i in 0...k do
        a = 0
        while a == 0
             a = rand(n)
        end
        if (!miller_rabin_pass(a, n))
            return false
        end
    end
    return true
end

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

Finally, some test code:

<<MillerRabin.rb>>=
Modular Exponent
Miller-Rabin pass
Miller-Rabin
if ARGV[0] == "test"
    n = ARGV[1].to_i
    puts (miller_rabin(n) ? "PRIME" : "COMPOSITE")
elsif ARGV[0] == "genprime"
    nbits = ARGV[1].to_i
    while true
        p = rand(2**nbits)
        if (p % 2 == 0)
            next
        elsif (p % 3 == 0)
            next
        elsif (p % 5 == 0)
            next
        elsif (p % 7 == 0)
            next
        end
        if miller_rabin(p)
            puts "#{p}"
            break
        end
    end
end

We've implemented here the modular_exponent.

<<Modular Exponent>>=
def modPow(x, r, m)
    y = r
    z = x
    v = 1
    while y > 0
        u = y % 2
        t = y / 2
        if u == 1
            v = (v * z) % m
        end
        z = z * z % m
        y = t
    end
    return v
end

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)
    next
elsif (p % 3 == 0)
    next
elsif (p % 5 == 0)
    next
elsif (p % 7 == 0)
    next
end

Here's some sample output:

$ ruby MillerRabin.rb test 516119616549881
PRIME
$ ruby MillerRabin.rb test 516119616549887
COMPOSITE
$ ruby MillerRabin.rb genprime 128
180226120047981445961500647274654519069
$ ruby MillerRabin.rb genprime 512
13294493160497232754943474450566763205781118232050806057882567789891631497772193\
268117571250170820851303260299628048340680597899360189633321043490933876713
$ ruby MillerRabin.rb genprime 1500
10607946794507398509197581673861929691947409838079439043722366281505052012659598\
20432361317748417405295826069160006004851781833481144992819835391487930979232076\
38035628804533994614946399143267306714617580824605075791836092123466680204438781\
16101151948165140759277188816046069504227443989047303572737109011537143396142206\
09172997967396676792816755270580692276847371160448864889127088914352989232555069\
0183514586128482948520257065750805597799560661164431

Download code