Miller-Rabin primality test

Remco Bloemen

2014-02-27, last updated 2014-03-04

As an exercise I implemented the Miller-Rabin primality test in plain C++. It turns out this algorithms lends itself for a true festival of operators, so I couldn’t resist making the code very dense:

bool MillerRabin(uint64 n, uint64 k)
{
    if(n == k) return true;
    uint64 s, d, b, e, x;

    // Factor n-1 as d 2^s
    for(s = 0, d = n - 1; !(d & 1); s++)
        d >>= 1;

    // x = k^d mod n using exponentiation by squaring
    // The squaring overflows for n >= 2^32
    for(x = 1, b = k % n, e = d; e; e >>= 1) {
        if(e & 1)
            x = (x * b) % n;
        b = (b * b) % n;
    }

    // Verify k^(d 2^[0…s-1]) mod n != 1
    if(x == 1 || x == n-1)
        return true;
    while(s-- > 1) {
        x = (x * x) % n;
        if(x == 1)
            return false;
        if(x == n-1)
            return true;
    }
    return false;
}

That’s it! Now, the algorithm will overflow for n232n \geq 2^{32} but I intend to use it to find 64-bit primes, so I’ll have to work on that.

Interestingly, the “primes page” claims that (source):

If n < 4,759,123,141 is a 2, 7 and 61-SPRP, then n is prime.

(A number n is called “kk-SPRP” or “strong probable-prime base kk” if it passes Miller Rabin with this k)

Since 4 759 123 141 is greater that 2322^{32} we can use this to implement a fast and exact function to determine whether some 32 bit number is prime. I special cased the first several primes to improve performance and made the whole thing into an incomprehensible operator soup, just like the Miller-Rabin implementation:

bool isPrime(uint32 n)
{
    return (n>73&&!(n%2&&n%3&&n%5&&n%7&&
    n%11&&n%13&&n%17&&n%19&&n%23&&n%29&&
    n%31&&n%37&&n%41&&n%43&&n%47&&n%53&&
    n%59&&n%61&&n%67&&n%71&&n%73))?false:
    MillerRabin(n, 2)&&MillerRabin(n, 7)
    &&MillerRabin(n, 61);
}

64-bits

In the previous post I showed an implementation of Miller-Rabin that overflowed when the number was more that 32 bits. In this post I present a way to fix this.

I could remember an instruction I used before, mulq which does 128 bit multiplication. It multiplies two 64 bit number and stores the least significant 64 bits in rax and the significant half in rdx.

Now we still need to reduce this number modulo m. Originally I planned to create a table of 2nmodm2^n \mod m, use the bits of rdx to add them and then add rax % m. But after consulting this x86-64 manual I found out about divq, which takes a 128 bit number stored in rdx:rax, a 64 bit number mm and produces the quotient and remainder. The whole thing can be coded as:

This should work fast, but I still need to look into Montgomery reduction. Perhaps I should benchmark the two.

Modular exponentiation can now be implemented using exponentiation by squaring:

And the Miller-Rabin test becomes:

Together with the prime sieve I created earlier a fast probablilistic prime checking procedure can be developed:

This algorithm correctly identifies the first ten primes less than 2642^{64}. But then again, so does a single Miller-Rabin round with K=17K= 17. Perhaps I should implement the Baillie-PSW algorithm, there are rumours that it has no false positives for 64 bits!