deterministic variant of Miller Rabbin method

output.txt

primality test using primality_miller_rabin(int64_t)
L = 100000000 : R = 100100000
test 1 duration : 0.861 s
test 2 duration : 0.853 s
test 3 duration : 0.856 s
test 4 duration : 0.844 s
test 5 duration : 0.84 s
test 6 duration : 0.835 s
test 7 duration : 0.839 s
test 8 duration : 0.839 s
test 9 duration : 0.838 s
test 10 duration : 0.836 s
average of 10 test(s) : 0.8441 s

primality_miller_rabin.cpp

typedef int64_t LL;

LL mulMod(LL a, LL b, LL mod) {  //computes a * b % mod
    LL r = 0;
    a %= mod, b %= mod;
    while (b) {
        if (b & 1) r = (r + a) % mod;
        b >>= 1, a = (a << 1) % mod;
    }
    return r;
}

LL powMod(LL a, LL n, LL mod) {  //computes a^n % mod
    LL r = 1;
    while (n) {
        if (n & 1) r = mulMod(r, a, mod);
        n >>= 1, a = mulMod(a, a, mod);
    }
    return r;
}

// deterministic variant of Miller Rabbin method
bool primality_miller_rabin(LL n) {
    const int pn = 3, p[] = {2, 3, 5};
    /*
        if n < 2,047, it is enough to test p = 2
        if n < 1,373,653, it is enough to test p = 2 and 3
        if n < 9,080,191, it is enough to test p = 31 and 73
        if n < 25,326,001, it is enough to test p = 2, 3, and 5
        if n < 3,215,031,751, it is enough to test p = 2, 3, 5, and 7
        if n < 4,759,123,141, it is enough to test p = 2, 7, and 61
        if n < 1,122,004,669,633, it is enough to test p = 2, 13, 23, and 1662803
        if n < 2,152,302,898,747, it is enough to test p = 2, 3, 5, 7, and 11
        if n < 3,474,749,660,383, it is enough to test p = 2, 3, 5, 7, 11, and 13
        if n < 341,550,071,728,321, it is enough to test p = 2, 3, 5, 7, 11, 13, and 17
        if n < 3,825,123,056,546,413,051, it is enough to test p = 2, 3, 5, 7, 11, 13, 17, 19, and 23
        if n < 18,446,744,073,709,551,616 = 264, it is enough to test p = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37
        if n < 318,665,857,834,031,151,167,461, it is enough to test p = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37
        if n < 3,317,044,064,679,887,385,961,981, it is enough to test p = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, and 41
    */
    for (int i = 0; i < pn; ++i)
        if (n % p[i] == 0) return n == p[i];
    if (n < p[pn - 1]) return false;
    LL s = 0, t = n - 1;
    while (~t & 1)
        t >>= 1, ++s;
    for (int i = 0; i < pn; ++i) {
        LL pt = powMod(p[i], t, n);
        if (pt == 1) continue;
        bool ok = false;
        for (int j = 0; j < s && !ok; ++j) {
            if (pt == n - 1) ok = true;
            pt = mulMod(pt, pt, n);
        }
        if (!ok) return false;
    }
    return true;
}