Problem: Determine if a 32-bit number is prime (deterministically)

Solution: (in C++)

// Bases to test. Using the first 4 prime bases makes the test deterministic
// for all 32-bit integers. See https://oeis.org/A014233.
int64_t bases[] = {2, 3, 5, 7};

inline int countTrailingZeros(uint64_t n) {
  if (n == 0)
    return 64;
  return __builtin_ctzll(n);
}

int64_t modularExponentiation(int64_t base, int64_t exponent, int64_t modulus) {
  int64_t res = 1;
  int64_t b = base % modulus;
  int64_t e = exponent;

  while (e > 0) {
    if (e & 1) {
      // Doesn't overflow because we assume 32-bit integer inputs
      res = (res * b) % modulus;
    }
    b = (b * b) % modulus;
    e >>= 1;
  }
  return res;
}

bool isPrime(int64_t n) {
  if (n < 2) return false;
  if (n < 4) return true;
  if (!(n & 1)) return false;

  int64_t d = n - 1;
  unsigned s = countTrailingZeros(d);
  d = d >> s;

  for (uint64_t a : bases) {
    if (n <= a) break;
    int64_t x = modularExponentiation(a, d, n);
    if (x == 1 || x == n - 1)
      continue;
    bool composite = true;
    for (unsigned r = 1; r < s; ++r) {
      // Doesn't overflow because it is at most n < 32 bits
      x = (x * x) % n;
      if (x == n - 1) {
        composite = false;
        break;
      }
    }
    if (composite)
      return false;
  }
  return true;
}

Discussion: In the late 1980’s Gary Miller and Michael Rabin came up with their now-famous Miller-Rabin primality test. See Wikipedia and my 2013 Program Gallery entry. It was notable in that it provided a randomized algorithm to check if an integer is prime, which had a tunable parameter to increase the probability of being correct.

The intervening 40 years has seen a huge body of research improving both randomized and deterministic primality testing. In 2002, Agrawal, Kayal, and Saxena found a condition-free deterministic polynomial-time algorithm, and the Ballie-PSW test, designed around the same time as Miller-Rabin, is still often used today in conjunction with Miller-Rabin. One of my favorite papers showing the importance of doing primality testing properly is in cryptography, Albrecht et al’s 2018 paper, Prime and Prejudice: Primality Testing Under Adversarial Conditions.

In relation to my work on homomorphic encryption, I found myself needing to generate a list of 40 primes, all roughly 32-bit in size, with particular properties. I stumbled across OEIS A014233, through which I learned about the study of strong pseudoprimes.

A strong pseudoprime is a composite number that passes Miller-Rabin’s deterministic test for a particular prime base $a$. That is, a number $n$ of the form $d \cdot 2^s + 1$ such that $a^d \equiv 1 \mod n$ or $a^{d \cdot 2^r} \equiv -1 \mod n$ for some $0 \leq r

If you test multiple bases on the same input, you’ll hedge against hitting small strong pseudoprimes. Testing 2, 3, and 5, means the smallest pseudoprime that confounds your test is 25326001. The code above demonstrates that if you add 7 to this list, you get a deterministic test for all 32-bit integers.

To the best of my knowledge, the idea to track the growth rate of strong pseudoprimes for the purpose of fast primality testing was first published in a 1980 Mathematics of Computation journal paper by Carl Pomerance, J. L. Selfridge and Samuel S. Wagstaff, Jr., titled “The Pseudoprimes to $25 \cdot 10^9$.”

The SPRP bases website has a nice little competition: who can find the best set of bases (not necessarily prime) for deterministic primality testing? For each cardinality (number of bases), the site lists the set of bases that produces the largest minimal pseudoprime to those bases. For covering all 64-bit integers, Jim Sinclair found this set of 7 bases.

For 32-bits, it can be done with three 64-bit bases

though the code above would need to be modified to account for overflow.

Performance-wise, the naive implementation of deterministic Miller-Rabin is decently fast. It can test primality of all 32-bit numbers in about 2 minutes on a Macbook (single-threaded). That said, there are far faster implementations, such as Kim Walisch’s primesieve, which can generate all 32-bit primes in 60 milliseconds. Notably, it does not use deterministic Miller-Rabin, opting instead for a sieve-based method with lots of cache optimizations and multithreading.

I am not a CPU performance tuning expert, but if you are it might be fun to try to beat that runtime using this method. The SPRP bases also includes a section where they discuss using hashing to reduce the compute cost of the compositeness test in Miller-Rabin, which seems like it would be a useful component. I asked Kim Walisch, the author of primesieve, and they replied that they had not tried a deterministic Miller-Rabin variant.

The code above is also on GitHub.