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