# Miller Rabin primality test (for 32 bit)

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);
}