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 kn

Complexity: O(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 apa(modp)

or alternatively

ap11(modp)

Proofs of this theorem can be found here

Some examples

351811(mod5)3111590491(mod11)

The converse of this theorem is not always true

If an11(modn) for some value of 0<a<n then n is prime

an example:

556111(mod561) but 561=31117 

but:

35611375(mod561)

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 ap11(modp) 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 a1,a2,,ai 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 ap12±1(modp)

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

a(2q+1)11(modp)

which means that

a2q10(modp)

this can be factored as

(aq1)(aq+1)0(modp)

therefore aq 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=(p1)2

Expressing Euler’s primality test formally:

If a(n1)/2±1(modn) where gcd(a,n)=1 then n must be a composite number for one of the following reasons:

  • if an11(modn) then n must be composite by Fermat’s Little Theorem
  • if an11(modn) then n must be composite because a(n1)/2 which is a square root of an1(modn) must fulfill the following equivalence a(n1)/2±1(modn) which is a condradiction to the statement above

This test also has some false positives e.g.

3(3411)/21(mod341) but 341=1131 

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 an1 it looks at the sequence of square roots/powers of two derived from an1

Let 2s be the largest power of 2 that divides n1, then n1=2sq for some odd integer q, the sequence of powers of two that divide n1 is

20,21,,2iwhere 0is

We know from Euler’s primality test that if an11(modn) then a(n1)/2±1(modn), let’s say that a(n1)/21(modn) then also because of Euler’s primality test a(n1)/22±1(modn), what this says is that as long as we can take the square root of some a(n1)/2i1(modn) the result must be ±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 an12i(modn) i.e. when n12i is no longer divisible by 2 which is exactly the number q, for this base case we’re sure of something, if aq±1(modn) then it means that it’s the square root of a2q1(modn) (obviously 2qn1 because n1 is even and must be divisible by at least 2)

If aq±1(modn) we have to analyze a2q(modn) and there’re three possible outcomes:

  • a2q1(modn) which by Euler’s primality test implies that aq±1(modn) which contradicts the statement above, therefore n is composite
  • a2q1(modn) which by Euler’s primality test implies that it’s the square root of some a2iq (where 0<i<s1), which will eventually become an11(modn) by successive squaring, therefore we can say that n is a probable prime
  • a2q±1(modn) 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;
}