100% deterministic Miller-Rabin primality test for N < 2^32

is_prime_u32.cpp

#include <iostream>
#include <cstdint>
#include <gmpxx.h>
#include <cassert>
#include <omp.h>

// Single-round Miller-Rabin test using GMP modular exponentiation
inline bool miller_rabin_test(uint32_t n, uint32_t a, uint32_t d, int s) {
    mpz_class base(a), exp(d), mod(n), x;
    
    // x = (a^d) % n
    mpz_powm(x.get_mpz_t(), base.get_mpz_t(), exp.get_mpz_t(), mod.get_mpz_t());

    if (x == 1 || x == n - 1) return true;

    for (int r = 1; r < s; r++) {
        // x = (x^2) % n
        mpz_powm_ui(x.get_mpz_t(), x.get_mpz_t(), 2, mod.get_mpz_t());
        if (x == n - 1) return true;
    }
    return false;
}

// 100% Deterministic Miller-Rabin Primality Test for N < 2^32
inline bool is_prime_u32(uint32_t n) {
    if (n < 2) return false;
    if (n == 2 || n == 3 || n == 5 || n == 7) return true;
    if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0) return false;

    uint32_t d = n - 1;
    int s = 0;
    while (d % 2 == 0) {
        d /= 2;
        s++;
    }

    // Proven base set for n < 2^32 (Jaeschke, 1993)
    static const uint32_t bases[] = {2, 7, 61};

    for (uint32_t a : bases) {
        if (n <= a) break;
        if (!miller_rabin_test(n, a, d, s)) return false;
    }

    return true;
}

int main() {
    const uint64_t N = 1ULL << 32;
    char *p = static_cast<char*>(malloc(N));
    assert(p);

    std::cout << "Initializing sieve..." << std::endl;
    for (uint64_t i = 0; i < N; ++i) p[i] = 1;
    p[0] = p[1] = 0;

    for (uint64_t i = 2; i * i < N; ++i) {
        if (p[i]) {
            for (uint64_t j = i * i; j < N; j += i) {
                p[j] = 0;
            }
        }
    }

    double start_time = omp_get_wtime();

    #pragma omp parallel
    {
        // Print thread count once
        #pragma omp single
        {
            std::cout << "Starting OpenMP verification across "
                      << omp_get_num_threads() << " threads..." << std::endl;
        }

        // Schedule belongs here on omp for:
        #pragma omp for schedule(guided)
        for (uint64_t i = 0; i < N; ++i) {
            bool expected = (p[i] == 1);
            bool actual = is_prime_u32(static_cast<uint32_t>(i));

            if (actual != expected) {
                #pragma omp critical
                {
                    std::cerr << "Mismatch found at index: " << i << std::endl;
                }
                assert(actual == expected);
            }
        }
    }

    double end_time = omp_get_wtime();

    std::cout << "All 2^32 tests passed successfully in " 
              << (end_time - start_time) << " seconds!" << std::endl;

    free(p);
    return 0;
}

Developed in longer chat with Gemini, includes proof of correctness.

$ g++ -O3 -Wall -Wextra -pedantic is_prime_u32.cpp -lgmpxx -lgmp -fopenmp
$ nproc
32
$ OMP_NUM_THREADS=16 ./a.out
Initializing sieve...
Starting OpenMP verification across 16 threads...
All 2^32 tests passed successfully in 21.8784 seconds!
$ ./a.out
Initializing sieve...
Starting OpenMP verification across 32 threads...
All 2^32 tests passed successfully in 16.2224 seconds!
$ 

$ grep "model name" /proc/cpuinfo | uniq -c
     32 model name	: AMD Ryzen 9 9950X 16-Core Processor
$ 
$ echo 0 | sudo tee /proc/sys/kernel/perf_event_paranoid
[sudo] password for hermann: 
0
$ OMP_NUM_THREADS=16 perf stat ./a.out 2>&1 | egrep "(elapsed|GHz)"
 1.757.434.277.570      cycles                           #    4,858 GHz                       
      38,321276381 seconds time elapsed
$ 
$ perf stat ./a.out 2>&1 | egrep "(elapsed|GHz)"
 2.339.569.609.756      cycles                           #    4,389 GHz                       
      32,392859350 seconds time elapsed
$ 

Used to build a prime list of biggest primes $&lt;2^{32}$ here:
https://github.com/Hermann-SW/uni-heidelberg/blob/main/scripts/proth_prover.cpp#L205-L218

    // 1. Build the prime list downwards from prime 4294967291
    double target_log = 2.0 * mpz_sizeinbase(N.get_mpz_t(), 2) * log(2.0);
    std::vector<uint64_t> primes_list;
    double current_log = 0.0;

    uint32_t p = UINT_MAX;
    while (current_log < target_log) {
        while (!is_prime_u32(p)) {
            p -= 2;
        }
        primes_list.push_back(p);
        current_log += log(p);
        p -= 2;
    }

Doing that is useful to find a list of primes whose product is greater than N^2 for very big number N (1000s of decimal digits). N % primes_list[i] are uint32_t by construction, so a square of that fits into uint64_t, and after mod the number is uint32_t again. Doing that from highest uint32_t downwards instead the normal 2, 3, 5, ... has the advantage that minimal number of uint32_t primes are needed for a given very big N.

So how far does this restriction to uint32_t primes help?
There are 203,280,221 primes up to $2^{32}$:

? primepi(2^32-1)
203280221
? 

PARI/GP has builtin primorial function, similar to factorial but only multiplying primes:

? [7!, factor(7!), 7#, factor(7#)]
[5040, [2, 4; 3, 2; 5, 1; 7, 1], 210, [2, 1; 3, 1; 5, 1; 7, 1]]
? 

So I wanted to determine the product of all above determined 203,280,221 primes, better the number of decimal digits of (needs 25GB resident memory to compute):

? log((2^32-1)#)/log(10)
1865246805.5323464684394081391321348541
?

So 1,865,246,805 is more than a billion digits for $N^2$, and number up to N with up to 932,623,402 decimal digits can be processed by that Residual Number System (RNS). Definitely good enough for 100 million decimal digit prime proof price money (https://www.eff.org/awards/coop).

181,918,401 primes starting with $2,3,5,\dots$ are needed for 1,660,000,000 decimal digits, leaving more than 200,000,000 decimal digits for $N^2$ or 100,000,000 decimal digits for N:

? T=166*10^7*log(10);
? S=0;forprime(p=2,oo,S+=log(p);if(S>T,print(primepi(p)," ",S);break));
181918401 3822291258.1464946014329212775077043265
? 
? S/log(10)
1660000001.6400604596343666596203056944
? 

Now the difference halved confirms that the memaining high primes product has enough digits:

? (1865246805-1660000001)\2
102623402
? 

And how much primes are needed for the product?
"Only" 21,361,820 primes:

? 203280221-181918401
21361820
? 

All calculations done for CPU with unit64_t/uint32_t.

Sounds that such RNS computations should be done on GPUs with thousands of cores, or even on many. In basement I have ten AMD gfx906 chip GPUs (8x Instinct MI50, 1x Radeon Pro VII, 1x Radeon VII) with 3,840 cores each, so 38,400 in total ...

Those GPUs do not have INT64, but FP64 (and provide >58 TFLOPS FP64 in total). To keep integers represented without precision loss, for that maximal prime precprime(sqrt(2^53)) would need to be used, again downwards as before.

? precprime(sqrt(2^53))
94906249
? 

That restriction leads to smaller maximal numbers to process, only 41.2 million decimal digits:

? S=0;forprime(p=2,94906249,S+=log(p));S/log(10)
41212381.378096148338898172464880684168
? 

Just a bit more than currently known biggest prime number $M_{52}$ with 41,024,320 decimal digits.

This reduced maximal digit number has one big advantage — "only" 5,484,598 different primes instead of 21.3 million:

? primepi(94906249)
5484598
?