Uses the deterministic Miller test and can test primality for values up to 2^32.
CountPrimes will count primes up to and including a value.
CountPrimesParallel will count all primes up to an exponent of 2 in parallel.
#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);
}
Top comments (0)