소수 생성 파트 2: 소수성을 위한 결정론적 밀러 테스트

#include <iostream>
#include <ppl.h>
#include <iomanip>
#include <Windows.h>

//Primes to test against
int witness[] = { 2, 7, 61 };

//number of primes up to exponents of 2
int pp[] = { 1, 2, 4, 6, 11, 18, 31, 54, 97, 172, 309, 564, 1028, 1900, 3512, 6542, 12251, 23000, 43390, 82025, 155611, 295947, 564163, 1077871, 2063689, 3957809, 7603553, 14630843, 28192750, 54400028, 105097565, 203280221 };

uint64_t Powmod(uint64_t n, uint64_t e, uint64_t m)
{
    uint64_t x = 1;
    while (e > 0)
    {
        if (e & 1)
            x = (x * n) % m;

        n = (n * n) % m;
        e >>= 1;
    }
    return x;
}

// See https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test.
// A direct translation of the deterministic Miller test, NOT the probabilistic variant.
bool MillerTest(uint64_t n)
{
    uint64_t d = n - 1, r = 0;

    while (!(d & 1))
    {
        r++;
        d >>= 1;
    }

    for (int i = 0; i < 3; i++)
    {
        uint64_t a = witness[i];
        uint64_t x = Powmod(a, d, n);

        //triggered if prime is in witness list
        if (x == 0)
            return true;

        if (x == 1 || x == (n - 1))
            continue;

        uint64_t j;

        for (j = 1; j < r; j++)
        {
            x = (x * x) % n;

            if (x == (n - 1))
                break;
        }

        if (j == r)
            return false;
    }
    return true;
}

void CountPrimes(uint64_t limit)
{
    uint64_t count = 0;
    time_t t1 = time(0);

    for (uint64_t i = 3; i <= limit; i += 2)
        if (MillerTest(i))
            count++;

    std::cout << "There are " << count + 1 << " primes up to " << limit << " " << time(0) - t1 << " seconds." << std::endl;


}

void CountPrimesParallel(uint64_t exponenetlimit)
{
    for (uint64_t i = 3; i <= exponenetlimit; i++)
    {
        uint64_t count = 0;
        time_t t1 = time(0);

        const uint64_t x = exp2l(i) / 2;
        std::atomic_uint count2 = 0;

        concurrency::parallel_for(uint64_t(1), x, [&](uint64_t a) {
            if (MillerTest(a * 2 + 1))
                count2++;
            });

        std::cout << "Exponent " << i << " " << x * 2 << " Result " << count2 + 1 << " Expected " << pp[i - 1] << " " << time(0) - t1 << " seconds." << std::endl;
    }
}


int main()
{
    CountPrimes(512);

    //2^32 is max value that can be tested
    CountPrimesParallel(32);
}