Primality Test

A prime number is a natural number greater than $1$ which has no positive divisors other than $1$ and itself

Naive test

Let $n$ be the number we want to check if is prime, if we find a natural number greater than $1$ that is a divisor of $n$ then $n$ is not a prime

• if a number $n$ is divisible by $k$ then $k \leq \sqrt{n}$

Complexity: $O(\sqrt{n})$

bool is_prime(int n) {
  if (n == 2) {
    // 2 is a prime number
    return true;
  }
  if (n == 1 || (n % 2 == 0)) {
    // 1 or any multiple of 2 is not a prime number
    return false;
  }
  for (int i = 3; i * i <= n; i += 2) {
    // check for any odd number < sqrt(n) if they are multiples of `n`
    if (n % i == 0) {
      return false;
    }
  }
  return true;
}

Erathostenes Sieve

If we have to make constants queries to check for numbers that are prime less than some number $n$ we can preprocess them using the Erathostenes Sieve and answer each query in $O(1)$

Fermat primality test

Fermat’s little theorem

If $a$ is an integer, $p$ a prime number where $0 < a < p$ then $$ a^p \equiv a \pmod{p} $$

or alternatively

$$ a^{p-1} \equiv 1 \pmod{p} $$

Proofs of this theorem can be found here

Some examples

The converse of this theorem is not always true

If $$ a^{n - 1} \equiv 1 \pmod{n} $$ for some value of $0 < a < n$ then $n$ is prime

an example:

but:

we can’t use the theorem directly to test if a number is prime since there’s a chance that the input is one of these special numbers (called Carmichael numbers) and the algorithm will give false positives e.g. $a = 5, p = 561$

what we can do is run the algorithm multiple times increasing the probability of finding a number $a$ such that $a^{p - 1} \not\equiv 1 \pmod{p}$ thus proving that $p$ is composite

// C++11
#include <random>

bool is_probably_prime(unsigned long long p, int iterations) {
  if (p == 2) {
    return true;
  }
  if (p % 2 == 0 || p == 1) {
    return false;
  }

  std::random_device rd;
  std::mt19937 engine(rd());
  std::uniform_int_distribution<long long> dis(2, p - 2);
  while (iterations--) {
    // choose an integer between 2 and n-2
    long long a = dis(engine);
    if (binary_exponentiation_modulo_m(a, p - 1, p) != 1) {
      return false;
    }
  }
  return true;
}

No matter how many iterations we use in the algorithm above there’s a chance that for each $a_1, a_2, \ldots, a_i$ Fermat’s little theorem holds true even though that the input is composite therefore this test is not used in practice

Euler primality test

Euler primality test is an improvement over the Fermat primality test because it adds another equality condition that a prime number must fulfill, assuming that $p$ is a prime number and $a$ is an integer where $0 < a < p$ then

If $a$ is an integer, $p$ a prime number where $0 < a < p$, $p > 2$ then $$ a^{\tfrac{p - 1}{2}} \equiv \pm 1 \pmod{p} $$

The motivation to this definition comes to the fact that any prime $> 2$ is an odd number, then the prime number can be expressed as $2q + 1$ where $q$ is an integer thus

which means that

this can be factored as

therefore $a^q$ is congruent to two possible values $1$ and $-1$. Going back to the definition of $q$, $2q + 1 = p$ we can find the value of $q$ as $q = \tfrac{(p - 1)}{2}$

Expressing Euler’s primality test formally:

If $a^{(n - 1) / 2} \not\equiv \pm 1 \pmod n$ where $gcd(a, n) = 1$ then $n$ must be a composite number for one of the following reasons:

• if $a^{n - 1} \not\equiv 1 \pmod{n}$ then $n$ must be composite by Fermat’s Little Theorem
• if $a^{n - 1} \equiv 1 \pmod{n}$ then $n$ must be composite because $a^{(n - 1) / 2}$ which is a square root of $a^{n - 1} \pmod{n}$ must fulfill the following equivalence $a^{(n - 1) / 2} \equiv \pm 1 \pmod n$ which is a condradiction to the statement above

This test also has some false positives e.g.

Miller-Rabin primality test

The Miller-Rabin primality test is quite similar to Euler’s primality test, but instead of looking at the square root of $a^{n - 1}$ it looks at the sequence of square roots/powers of two derived from $a^{n - 1}$

Let $2^s$ be the largest power of $2$ that divides $n - 1$, then $n - 1 = 2^s \cdot q$ for some odd integer $q$, the sequence of powers of two that divide $n - 1$ is

We know from Euler’s primality test that if $a^{n - 1} \equiv 1 \pmod{n}$ then $a^{(n - 1) / 2} \equiv \pm 1 \pmod{n}$, let’s say that $a^{(n - 1) / 2} \equiv 1 \pmod{n}$ then also because of Euler’s primality test $a^{(n - 1) / 2^2} \equiv \pm 1 \pmod{n}$, what this says is that as long as we can take the square root of some $a^{(n - 1) / 2^i} \equiv 1 \pmod{n}$ the result must be $\pm 1$ otherwise it’s a composite number by Euler’s primality test

The base case occurs when we cannot take the square root of some $a^{\tfrac{n - 1}{2^i}} \pmod{n}$ i.e. when $\tfrac{n - 1}{2^i}$ is no longer divisible by $2$ which is exactly the number $q$, for this base case we’re sure of something, if $a^q \equiv \pm 1 \pmod{n}$ then it means that it’s the square root of $a^{2q} \equiv 1 \pmod{n}$ (obviously $2q \leq n - 1$ because $n - 1$ is even and must be divisible by at least $2$)

If $a^q \not\equiv \pm 1 \pmod{n}$ we have to analyze $a^2q \pmod{n}$ and there’re three possible outcomes:

• $a^2q \equiv 1 \pmod{n}$ which by Euler’s primality test implies that $a^q \equiv \pm 1 \pmod{n}$ which contradicts the statement above, therefore $n$ is composite
• $a^2q \equiv -1 \pmod{n}$ which by Euler’s primality test implies that it’s the square root of some $a^{2^iq}$ (where $0 < i < s-1$), which will eventually become $a^{n - 1} \equiv 1 \pmod{n}$ by successive squaring, therefore we can say that $n$ is a probable prime
• $a^2q \not\equiv \pm 1 \pmod{n}$ which is the same as the statement above (therefore we have to keep analyzing the next element in the sequence)

// C++11
#include <random>

bool miller_rabin_primality_test(long long a, long long n) {
  int s = 0;
  long long q = n - 1;
  while (q % 2 == 0) {
    q /= 2;
    s += 1;
  }
  long long m = binary_exponentiation_modulo_m(a, q, n);
  if (m == 1 || m == n - 1) {
    // base case a^q ≡ 1 (mod n)
    return true;
  }
  for (int i = 0; i < s; i += 1) {
    // a^{2^iq} (mod n)
    m = (m * m) % n;
    if (m == n - 1) {
      return true;
    }
  }
  return false;
}

bool is_probably_prime(long long p, int iterations) {
  // NOTE: test of the primes 2 and 3 because of
  // the distribution limits (p, p - 2)
  if (p == 2 || p == 3) {
    return true;
  }
  if (p % 2 == 0 || p == 1) {
    return false;
  }
  std::random_device rd;
  std::mt19937 engine(rd());
  std::uniform_int_distribution<long long> dis(2, p - 2);
  while (iterations--) {
    // choose an integer between 2 and n-2
    long long a = dis(engine);
    if (!miller_rabin_primality_test(a, p)) {
      return false;
    }
  }
  return true;
}