Primality Test

This post introduces some popular algorithms for checking if a number is prime or not.

Deterministic algorithms

Simplest algorithm

A prime number is a natural number greater than 1 and has only two divisors: 1 and itself. Following this definition, the most basic algorithm for checking if n is prime would be:

Check if any number between 2 and n-1 is a divisor of n. If there exists one, n is not prime, otherwise it is.

Sample implementation in C++:

bool isPrime(int n) {
  for (int i = 2; i <= n - 1; ++i) {
    if (n % i == 0) {
      return false;
    }
  }
  return true;
}

Time complexity of this algorithm is O(n), not fast enough for n up to 1e9.

Better algorithm

It d is a divisor of n, \dfrac{n}{d} must also be a divisor of n. It’s easy to see that either d or \dfrac{n}{d} is smaller than \sqrt{n}. Therefore, we only need to check for numbers between 2 and \sqrt{n}.

Sample implementation in C++:

bool isPrime(int n) {
  for (int i = 2; i * i <= n; ++i) {
    if (n % i == 0) {
      return false;
    }
  }
  return true;
}

This algorithm has a complexity of \sqrt{n} and is often used in practice. It can check for n up to 1e14. There are some optimizations available (e.g. don’t check even numbers other than 2), but the complexity is essentially the same. I won’t discuss them here.

Probabilistic algorithms

Fermat primality test

Fermat’s little theorem states that for a prime number p and a integer a not divisible by p:

a^{p-1} \equiv 1 \bmod p

To test whether some number p is prime, we can pick random integers a not divisible by p and check if the equation holds.

If the equation does not hold, we know for sure that p is composite. If the equation holds for some value of a, p is called a probable prime.

In practice, a is usually chosen randomly from the interval [2, p-2]. The test will be repeated multiple times for different values of a.

Sample implementation in C++:

int powmod(int a, int b, int mod) {
  int res = 1;
  for (; b; (a *= a) %= mod, b >>= 1)
    if (b & 1) (res *= a) %= mod;
  return res;
}

bool fermat(int n, int iter = 5) {
  if (n < 4) return (n == 2 || n == 3);

  for (int i = 1; i <= iter; i++) {
    int a = 2 + rand() % (n - 3);
    if (powmod(a, n - 1, n) != 1)
      return false;
  }
  return true;
}

powmod(a, b, c) uses binary exponentiation to effectively calculate a^{b} \mod c.

There’s one big problem though: a Carmichael number will pass this test for every base a co-prime to the number, even though it is not actually prime. Due to this, you shouldn’t use this test in your code or you will probably get hacked.

Miller–Rabin primality test

Miller–Rabin primality test extends upon the idea of Fermat primality test.

Let p be a prime number larger than 2. p must be odd, so p-1 is even. Factor out all powers of 2 in p-1, we have:

p-1 = 2^{s} \cdot d, with d odd

We use it to factorize the equation of Fermat’s little theorem:

\begin{array}{rl} &a^{n-1} \equiv 1 \bmod p \\ \Leftrightarrow &a^{2^s d} – 1 \equiv 0 \bmod p \\ \Leftrightarrow &(a^{2^{s-1} d} + 1) (a^{2^{s-1} d} – 1) \equiv 0 \bmod p \\ \Leftrightarrow &(a^{2^{s-1} d} + 1) (a^{2^{s-2} d} + 1) (a^{2^{s-2} d} – 1) \equiv 0 \bmod p \\ \quad\vdots &\\ \Leftrightarrow &(a^{2^{s-1} d} + 1) (a^{2^{s-2} d} + 1) \cdots (a^{d} + 1) (a^{d} – 1) \equiv 0 \bmod p \\ \end{array}

If p is prime, p must divide one of these factors. That’s exactly what we check in the Miller–Rabin primality test. More specifically, given a number n and a base 2 \le a \le n-2, we check if either:

  • a^d \equiv 1 \bmod n holds or,
  • a^{2^r d} \equiv -1 \bmod n holds for some 0 \le r \le s-1

If we found a base a that doesn’t satisfy any of these conditions, again we know for sure that p is composite. Otherwise, n is called a strong probable prime. If n is indeed composite, a is called a strong liar.

It can be proved that at most 25\% of bases a can be strong liars. We repeat the test multiple times with different bases. This makes the probability that a composite number gets declared as a strong probable prime very small. Besides, there’s nothing like the Carmichael numbers.

Overall, the Miller–Rabin test has much better accuracy than the Fermat test and can be safely used in competitive programming.

Deterministic version

It’s possible to make the test deterministic by trying all possible bases a \le 2 \ln(n)^2. Still, this limit is pretty large, so people have invested a lot of computational power in finding better limits.

For 32-bit integers, it suffices to check only the first 4 prime bases: 2, 3, 5, 7.

For 64-bit integers, it suffices to check only the first 12 prime bases: 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37.

Sample implementation in C++ for 64-bit integers:

using ull = unsigned long long;

// for 64 bit integer it is enough to check the first 12 prime bases
const ull BASE[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

ull mulmod(ull a, ull b, ull mod) {
  ull res = 0;
  for (; b; (a += a) %= mod, b >>= 1)
    if (b & 1) (res += a) %= mod;
  return res;
}

ull powmod(ull a, ull b, ull mod) {
  ull res = 1;
  for (; b; a = mulmod(a, a, mod), b >>= 1)
    if (b & 1) res = mulmod(res, a, mod);
  return res;
}

bool check_prime(ull n, ull a, ull d, int s) {
  ull x = powmod(a, d, n);
  if (x == 1 || x == n - 1) return true;
  for (int r = 2; r <= s; ++r) {
    x = mulmod(x, x, n);
    if (x == n - 1) return true;
  }
  return false;
}

bool miller_rabin(ull n) {
  if (n < 4) return (n == 2 || n == 3);
  ull d = n - 1;
  int s = 0;
  while ((d & 1) == 0) ++s, d >>= 1;
  for (auto a: BASE) {
    if (n == a) return true;
    if (!check_prime(n, a, d, s)) return false;
  }
  return true;
}
5 1 vote
Article Rating
Subscribe
Notify of
guest
0 Comments
Inline Feedbacks
View all comments
0
Would love your thoughts, please comment.x
()
x