Algoteka
Log in to submit your own samples!

Miller–Rabin Primality Test

Problem by oml1111
# Tech tags Title Creator Created date
1 0
2022-11-26 19:19
View all samples for this language (1 verified and 0 unverified)

Simple 64-bit implementation | C++ |

By oml1111 |
0 likes

A simple implementation of the Miller-Rabin primality test algorithm for 64 bit integers. Runs in \(O(k \log^2 n)\) thanks to the fact that modular addition is constant for 64-bit integers (\(k\) is the number of tests, \(n\) is the integer value).

Caveat: The value of \(2n\) must not overflow a 64-bit integer. For larger integers, you need to use some big-integer class with addition and modulo operations implemented

Code

#include<random>


long long modmult(long long a, long long b, long long mod) {  // Compute a*b % mod
    long long result = 0;
    while(b > 0) {
        if(b%2 == 1) {
            result = (result+a) % mod;
        }
        a = (a+a) % mod;
        b /= 2;
    }
    return result;
}


long long modpow(long long a, long long b, long long mod) {  // Compute a^b % mod
    long long result = 1;
    while(b > 0) {
        if(b%2 == 1) {
            result = modmult(result, a, mod);
        }
        a = modmult(a, a, mod);
        b /= 2;
    }
    return result;
}


// Miller–Rabin primality test
bool prime_test(long long n, int num_tests = 20) {  // Assumes n is an odd integer larger than 3
    static std::mt19937_64 randgen(0);  // Change seed if needed
    
    long long d = n-1;
    long long s = 0;
    while(d%2 == 0) {
        s++;
        d /= 2;
    }
    
    for(int test = 0; test < num_tests; test++) {
        long long a = randgen()%(n-3) + 2;  // Between 2 and v-2
        long long x = modpow(a, d, n);
        
        for(int i=0; i<s; i++) {
            long long y = modmult(x, x, n);
            if(y == 1 && x != 1 && x != n-1) {  // Nontrivial square root of 1 modulo n
                return false;  // (x+1)(x-1) divisible by n, meaning gcd(x+1, n) is a factor of n, negating primality
            }
            x = y;
        }
        if(x != 1) {
            return false;
        }
    }
    return true;  // Number is prime with likelihood of (1/4)^num_tests
}


int main() {
    prime_test(5127821565631733);
}

References

objects
std::mt19937_64 en.cppreference.com cplusplus.com
functions
std::mersenne_twister_engine::mersenne_twister_engine en.cppreference.com cplusplus.com
std::mersenne_twister_engine::operator() en.cppreference.com cplusplus.com

Problem Description

Implement the Miller-Rabin primality test algorithm. Use your implementation to test the primality of some number.

Miller–Rabin primality test - wikipedia.org

View sample discussion (0 comments)
View problem discussion (0 comments)