miller_rabin_probabilistic_test.cpp

int tt = 10, it = 1;
double total = 0.;
cout << "primality test using primality_miller_rabin_probabilistic(int64_t)" << endl;
cout << "Miller Rabin iteration count = " << it << endl;
cout << "L : " << L << " R : " << R << endl;
for(int k = 0, c = 0, t = clock(); k < tt; k++, c = 0, t = clock()) {
    for(int i = L; i <= R; i++)
        c += primality_miller_rabin_probabilistic(i, it);
    double time = (clock() - t) / (CLOCKS_PER_SEC + 0.);
    total += time;
    cout << "test " << k << " result : number of primes = " << c << " : time taken = " << time << endl;
}
cout << "average time for " << tt << " test(s) = " << total / tt << endl;

output_100000000_100100000.txt

primality test using primality_miller_rabin_probabilistic(int64_t)
Miller Rabin iteration count = 1
L : 100000000 R : 100100000
test 0 result : number of primes = 5411 : time taken = 0.874
test 1 result : number of primes = 5411 : time taken = 0.828
test 2 result : number of primes = 5411 : time taken = 0.833
test 3 result : number of primes = 5411 : time taken = 0.821
test 4 result : number of primes = 5411 : time taken = 0.819
test 5 result : number of primes = 5411 : time taken = 0.824
test 6 result : number of primes = 5411 : time taken = 0.812
test 7 result : number of primes = 5411 : time taken = 0.84
test 8 result : number of primes = 5411 : time taken = 0.816
test 9 result : number of primes = 5411 : time taken = 0.848
average time for 10 test(s) = 0.8315

output_10000000_10100000.txt

primality test using primality_miller_rabin_probabilistic(int64_t)
Miller Rabin iteration count = 1
L : 10000000 R : 10100000
test 0 result : number of primes = 6243 : time taken = 0.62
test 1 result : number of primes = 6242 : time taken = 0.6
test 2 result : number of primes = 6242 : time taken = 0.605
test 3 result : number of primes = 6241 : time taken = 0.596
test 4 result : number of primes = 6241 : time taken = 0.605
test 5 result : number of primes = 6243 : time taken = 0.588
test 6 result : number of primes = 6241 : time taken = 0.589
test 7 result : number of primes = 6241 : time taken = 0.588
test 8 result : number of primes = 6242 : time taken = 0.588
test 9 result : number of primes = 6242 : time taken = 0.593
average time for 10 test(s) = 0.5972

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

bool primality_miller_rabin_probabilistic(LL n, int k) {
    if(n % 2 == 0) return n == 2;
    if(n % 3 == 0) return n == 3;
    LL s = 0, d = n - 1;
    while(~d & 1) s++, d >>= 1;
    for(int i = 0; i < k; i++) {
        LL a = 2 + rand() % (n - 3);
        LL b = powMod(a, d, n);
        if(b != 1) {
            bool flag = true;
            for(int r = 0; r < s; r++, b = mulMod(b, b, n))
                if(b == n - 1) {
                    flag = false; break;
                }
            if(flag) return false;
        }
    }
    return true;
}

test.cpp

srand(time(NULL));

/* prime sieve starts */
LL U = 100000000;
double r = sqrt(U);
vector<bool> boolar(U, true);
for(int i = 3; i <= r; i += 2)
    if(boolar[i])
        for(int j = 3 * i, r = j - i; j < U; j += r)
            boolar[j] = false;
/* orime sieve ends */
vector<int> primepi(U, 1);
for(int i = 3; i < U; i++)
    if(i % 2 == 0)
        primepi[i] = primepi[i - 1];
    else
        primepi[i] = primepi[i - 1] + boolar[i];

cout << setprecision(15) << fixed;

int tt = 5;
for(int i = 1000000; i < 100000000; i += 1000000) {
    int L = i, R = L + 100000;
    int errcount = 0;
    double T = clock();
    for(int j = L; j <= R; j++) {
        bool ip = primepi[j] - primepi[j - 1];
        for(int k = 0; k < tt; k++)
            errcount += primality_miller_rabin_probabilistic(j, 2) & !ip;
    }
    cout << "{" << L << ", ";
    cout << (clock() - T) / (tt * CLOCKS_PER_SEC) << ", ";
    cout << errcount / (tt * (primepi[R] - primepi[L] + 0.)) <<"}, ";
}