# 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 $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 < 4,759,123,141 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 4 759 123 141 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);
}
```

# 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 $2^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 $m$ and produces the quotient and remainder.
The whole thing can be coded as:

```
/// Modular multiplication
/// @param a The first factor, a < m
/// @param a The second factor, b < m
/// @param m The modulus
/// @return The reduced product, a b mod m < m
uint64 mul(uint64 a, uint64 b, uint64 m)
{
// Perform 128 multiplication and division
uint64 q; // q = ⌊a b / m⌋
uint64 r; // r = a b mod m
asm("mulq %3;"
"divq %4;"
: "=a"(q), "=d"(r)
: "a"(a), "rm"(b), "rm"(m));
return r;
}
```

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:

```
/// Modular exponentiation
/// @param b The base, b < m
/// @param e The exponent
/// @param m The modulus
/// @returns The reduced power of a, a^b mod m
uint64 pow(uint64 b, uint64 e, uint64 m)
{
uint64 r = 1;
for(; e; e >>= 1) {
if(e & 1) r = mul(r, b, m);
b = mul(b, b, m);
}
return r;
}
```

And the Miller-Rabin test becomes:

```
/// Miller-Rabin probabilistic primality testing
/// @param n The number to test for primality
/// @param k The witness for primality
/// @returns True iff when n is a k-stong pseudoprime
bool MillerRabin(uint64 n, uint64 k)
{
// Factor n-1 as d*2^s
uint64 s = 0;
uint64 d = n - 1;
for(; !(d & 1); s++)
d >>= 1;
// Verify x = k^(d 2^i) mod n != 1
uint64 x = pow(k % n, d, n);
if(x == 1 || x == n-1) return true;
while(s-- > 1) {
// x = x^2 mod n
x = mul(x, x, n);
if(x == 1) return false;
if(x == n-1) return true;
}
return false;
}
```

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

```
const vector<uint64> small_primes = prime_sieve(1<<20);
/// Miller-Rabin probabilistic primality testing
/// @param n The number to test for primality
/// @returns False when n is not a prime
bool is_prime(uint64 n)
{
// Handle small primes fast
for(int i = 0; i < small_primes.size(); i++) {
uint64 p = small_primes[i];
if(n == p) return true;
if(n % p == 0) return false;
}
// Do a few Miller-Rabin rounds
for(int i = 0; i < 10; i++) {
if (!MillerRabin(n, small_primes[i])) return false;
}
return true;
}
```

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