# The Miller-Rabin primality test

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