% 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!