#include <iostream>
#include <ppl.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 };

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

//Count primes up to a specified unsigned 64 bit integer..
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;
}

//Count primes up to a specified power of 2 in parallel.
void CountPrimesParallel(uint64_t exponent)
{
    time_t t1 = time(0);
    const uint64_t x = exp2l(exponent) / 2;
    std::atomic_uint count = 0;

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

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

int main()
{
    //CountPrimes(512);

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