The Miller-Rabin primality test

% Remco Bloemen % 2014-05-07

The Miller-Rabin primality test

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

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

/// Miller-Rabin probabilistic primality testing
/// @param n The number to test for  primality
/// @param k The witness for primality
/// @returns True iff 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;
}

https://miller-rabin.appspot.com/

bool isPrime(uint64 n)
{
	// Small primes
	if(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)
		return false;

	// Jim Sinclair's Miller-Rabin base for 2^64
	if(!millerRabin(n, 2) || !millerRabin(n, 325) || !millerRabin(n, 9375))
	|| !millerRabin(n, 28178) || !millerRabin(n, 450775)
	|| !millerRabin(n, 9780504) || !millerRabin(n, 1795265022))
		return false;

	// Must be prime
	return true;
}