DEV Community

Antidisestablishmentarianism
Antidisestablishmentarianism

Posted on • Edited on

1 1

Prime Number Generation Part 2: Deterministic Miller Test For Primality

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);
}
Enter fullscreen mode Exit fullscreen mode

Qodo Takeover

Introducing Qodo Gen 1.0: Transform Your Workflow with Agentic AI

Rather than just generating snippets, our agents understand your entire project context, can make decisions, use tools, and carry out tasks autonomously.

Read full post →

Top comments (0)

Qodo Takeover

Introducing Qodo Gen 1.0: Transform Your Workflow with Agentic AI

Rather than just generating snippets, our agents understand your entire project context, can make decisions, use tools, and carry out tasks autonomously.

Read full post