Miller Rabin primality test (for 32 bit)

% Remco Bloemen % 2009-11-25, last updated 2014-03-03

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 n \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 < 4759123141 is a 2, 7 and 61-SPRP, then n is prime.

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

Since 4759123141 is greater that 2^{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);
}