Primality testing algorithms

In this post we're going to look at some algorithms for checking if a given number is either prime or composite.

Édouard Lucas, Leonhard Euler, Carl Pomerance and Pierre de Fermat

Prime numbers play an important role in public-key cryptography (e.g. RSA algorithm) which require fast generation of large random prime numbers that have around `600` digits (or more).

The simplest way to generate a prime number, is to pick a random odd integer and check its primality, repeating this process until a prime is found. This approach requires a fast algorithm for checking the primality of a number.

In 2002, Manindra Agrawal, Neeraj Kayal, and Nitin Saxena published the AKS primality test, unconditionally proving that the primality of a number can be checked in polynomial time, with an asymptotic time complexity of `Õ(log(n)^12)`.

Later on, in 2005, Carl Pomerance and Hendrik Lenstra improved the complexity of the AKS algorithm to run in just `Õ(log(n)^6)` steps.

While being of great theoretical importance, the AKS algorithm is not actually used in practice.

# Fermat primality test

For the most part of history, the fastest known algorithm for checking the primality of a number `p`, was by trying to divide it by smaller primes not exceeding `sqrt(p)`. If no prime divides it, then the number is prime. However, this is a very slow test for large `p`.

Everything started to change with Fermat's little theorem, which states that for any prime number `p` and any integer `a`:

`a^p ≡ a mod p`

If `gcd(a,p) = 1`, we can divide both sides by `a` and get the commonly used form:

`a^(p-1) ≡ 1 mod p`

This is a fast test due to the fact that modular exponentiation can be computed efficiently using the exponentiation by squaring algorithm.

However, it's not a perfect test, because, for every base `a`, there exist infinitely many composite numbers, known as Fermat pseudoprimes to base `a`, that also pass the test.

The list of Fermat pseudoprimes to base `2`, is given by A001567.

To make things worse, there also exist composite numbers `n`, behaving just like the primes, passing the Fermat test for every base `a` coprime to `n`. These numbers are known as Carmichael numbers (OEIS: A002997).

# Solovay-Strassen primality test

The Solovay-Strassen primality test is an improvement over Fermat's primality test, and is based on Euler's criterion, which states that for any prime number `p` and any integer `a`:

`a^((p-1)/2) ≡ (a/p) mod p`

where `(a/p)` is the Jacobi symbol.

By choosing several random bases `a` in the range `[2, n-1]`, we check if:

`a^((n-1)/2) ≡ (a/n) mod n`

If a given number `n` passes the test `k` times for random bases of `a`, then it's probably prime, with `2^(-k)` probability of actually being composite.

# Miller-Rabin primality test

The Miller-Rabin primality test is an improvement over the Solovay-Strassen primality test.

Let `n` be an odd prime number. By writing `n-1 = d * 2^s`, where `d` is odd, then for every base `a ∈ [2, n-1]`, either (1) or (2) is true:

`a^d ≡ 1 mod n`
`a^(d*2^r) ≡ -1 mod n`, for some `0 <= r < s`

The test is based on the contrapositiveness of the above statement. If we can find an integer `a ∈ [2, n-1]` such that :

`a^d ≢ 1 mod n`

and

`a^(d*2^r) ≢ -1 mod n`, for all `0 <= r < s`

then `n` is definitely composite.

By repeating the test `k` times, the probability of a composite being identified as prime is `4^(-k)`.

In 1976, Gary L. Miller noted that under the assumption of the Generalized Riemann Hypothesis (GRH) for Dirichlet characters, the test can be made deterministic by checking all the bases `a` in the range `2 <= a < C * log(n)^2`, for some constant `C`. In 1984, Eric Bach showed that we can use `C = 2`.

# Fibonacci primality test

The Fibonacci primality test uses the modular Fibonacci sequence `F_n`, where for any prime number `n`:

`F_(n - k) ≡ 0 mod n`

where `k = (5/n)` is the Jacobi symbol.

The list of odd Fibonacci pseudoprimes is given by A081264.

Research problem: find a composite number `n ≡ ±2 mod 5` that is both a Fibonacci pseudoprime and a Fermat pseudoprime to base `2`. Pomerance, Selfridge and Wagstaff offer `$620` for the first example or for a proof that no such number exists.

# Lucas primality test

The Lucas primality test is a generalization of the Fibonacci primality test, using the modular Lucas sequences `U_n(P,Q)` and `V_n(P,Q)`, where for any prime number `n`:

`V_n(P,Q) ≡ P mod n`

`U_(n-k)(P,Q) ≡ 0 mod n`

where `k = (D/n)` is the Jacobi symbol, for `D = P^2 - 4Q`, with `1 < P < n` and `Q < n` and `D != 0`.

The Lucas sequences can be computed efficiently, using the binary expansion of `n`.

For a fast algorithm, see the paper "Efficient computation of full Lucas sequences", by M. Joye and J.-J. Quisquater.

When `gcd(n, D) = 1`, a slightly faster algorithm can be used, described in the paper "On Lucas Sequences Computation", by Aleksey Koval.

# PSW primality test

The PSW primality test, named after Carl Pomerance, John Selfridge, and Samuel Wagstaff, uses two tests: a Miller-Rabin test to base `2` and a Lucas test to the `U_n(P,Q)` sequence, with a carefully chosen value of `D`.

Given an odd integer `n`, that is not a perfect power:

Perform a Miller-Rabin test to base `2`.
Find the first `P >= 1` for which the Jacobi symbol `((P^2+4)/n) = -1`.
If `U_(n+1)(P,-1) ≡ 0 mod n`, then `n` is probably prime.

No counter-examples are known to this test.

# Baillie-PSW primality test

The Baillie-PSW primality test, named after Robert Baillie, Carl Pomerance, John Selfridge, and Samuel Wagstaff, is a combination of two primality tests:

a Miller-Rabin test to base `2`
a strong Lucas test, with a carefully chosen value of `D`

The strong Lucas test can also be replaced with an extra-strong Lucas test.

In the strong Lucas test, the recommended value for `D` has the form:

`D = (-1)^k (2k + 1)`

with `(D/n) = -1`, for the smallest such integer `k >= 2`.

After finding `D`, set `P = 1` and `Q = (1-D)/4`, and express `n+1 = d * 2^s`, with `d` odd.

A strong Lucas pseudoprime must satisfy either (1) or (2):

`U_d(P,Q) ≡ 0 mod n`
`V_(d * 2^r)(P,Q) ≡ 0 mod n`, for some `0 <= r < s`

An extra-strong Lucas pseudoprime for a set of parameters `(P,Q)`, must satisfy either (1) or (2):

`U_d(P,Q) ≡ 0 mod n` and `V_d(P,Q) ≡ ±2 mod n`
`V_(d * 2^r)(P,Q) ≡ 0 mod n`, for some `0 <= r < s`

with `P >= 3` and `Q = 1` and `((P^2 - 4)/n) = -1`.

The lists of strong and extra-strong Lucas pseudoprimes are given by A217255 and A217719, respectively.

In the selection of parameters, if we encounter `(D/n) = 0`, where `D = P^2 - 4Q` and `1 < D < n`, then `n` is composite with `gcd(D,n)` being a non-trivial factor of `n`, allowing us to return early.

If `n` passes the Miller-Rabin test to base `2` and is also a strong Lucas pseudoprime, using the parameter selection described above, then `n` is almost surely a prime number. In fact, since its publication in `1980`, no counter-examples are known to this test.

# Quadratic Frobenius primality test

Another practical test, is the quadratic Frobenius primality test, developed by Jon Grantham, which also uses the modular Lucas sequences `U_n(P,Q)` and `V_n(P,Q)`, where for every prime number `n > 5`, the following congruences are satisfied:

`U_(n-k)(P, Q) ≡ 0 mod n`
`V_(n-k)(P,Q) ≡ 2 mod n`, if `k = 1`
`V_(n-k)(P,Q) ≡ 2Q mod n`, if `k = -1`

where `k = (D/n)`, with `D = P^2 - 4Q` and `gcd(n, 2 Q D) = 1` .

The test can be implemented efficiently, using the algorithm described by Aleksey Koval in the paper "On Lucas Sequences Computation".

In choosing the parameters `(P,Q)`, we can set `Q = 2` and search for the first odd integer `P >= 5`, satisfying `(D/n) = -1`.

No counter-examples are known to this test, using the selection of parameters described above.

# Tribonacci primality test

It's also possible to create primality testing methods of higher order, as illustrated by the following Tribonacci-like sequence:

`T(0) = 0`
`T(1) = 0`
`T(2) = 9`
`T(n) = T(n-1) + 3 T(n-2) + 9 T(n-3)`

where the first few terms are:

`T(n) = {0, 0, 9, 9, 36, 144, 333, 1089, 3384, 9648, 29601, 89001, 264636, 798048, ...}`

For prime numbers `n > 5`, we have:

`T(n) ≡ 0 mod n`, for prime `n > 5` congruent to `{1,3} mod 8`.
`T(n) ≡ 4 mod n`, for prime `n > 5` congruent to `{5,7} mod 8`.

This test can be implemented somewhat efficiently by using matrices and exponentiation by squaring:

`A = [[0, 3, 0],[0, 0, 3],[1,1,1]]`

If `n` is congruent to `{1,3} mod 8`, then:

`A^(n-1) ≡ [[1,0,0], [0,1,0], [0,0,1]] mod n`

If `n` is congruent to `{5,7} mod 8`, then:

`A^(n+1) ≡ [[4,2,3], [1,5,3], [1,2,6]] mod n`

When `n ≡ 1 mod 8`, the first few counter-examples to this test, are:

`88561, 107185, 162401, 221761, 226801, 334153, 410041, 665281, 825265, 1569457, 1615681, 2727649, ...`

When `n ≡ 3 mod 8`, several counter-examples include (not necessarily consecutive):

`80375707, 154287451, 267559627, 326266051, 478614067, 573183451, 643767931, 2433943891, 4297753027, ...`

When `n` is congruent to `{5,7} mod 8`, no counter-examples are known.

All known counter-examples are also Fermat pseudoprimes to base `2` and/or `3`.

A faster implementation of this method is possible by using quadratic integers, based on the following two congruences:

`(1 - sqrt(-2))^(n-1) ≡1 mod n`, if `n` is congruent to `{1,3} mod 8`

`(1 - sqrt(-2))^(n+1) ≡3 mod n`, if `n` is congruent to `{5,7} mod 8`

However, this method is more of a curiosity than a practical test.

# Pépin-Proth primality test

The Pépin test can be used for deterministically checking the primality of Fermat numbers `F_n = 2^(2^n)+1`:

`F_n` is prime if and only if `3^((F_n - 1)/2) ≡ -1 mod F_n`

The test itself is a variant of the Proth test, which is used for proving the primality of Proth numbers of the form `p = k * 2^n + 1` with `k` odd and `k < 2^n`, which is prime if and only if there exists an integer `a` such that:

`a^((p-1)/2) ≡ -1 mod p`

The Proth test can generalized to other numbers, by letting `n` to be a positive odd integer that is not a perfect power.

Algorithm:

find an integer `a` such that the Jacobi symbol `(a/n) = -1`.
if `n` is prime, then `a^((n-1)/2) ≡ -1 mod n`

However, there exist composite numbers `n` that also satisfy this property for some values of `a`.

To overcome this, we try to find `a` in the sequence `k, k+1, k+2, ...`, where `k` is a random integer less than `floor(sqrt(n))`. For better accuracy, we repeat the test several times. A random value of `k` is computed at each iteration.

If we find an integer `a` such that the Jacobi symbol `(a/n) = -1` and `a^((n-1)/2) ≢ -1 mod n`, then `n` is definitely composite.

# Lucas-Lehmer primality test

The Lucas-Lehmer test (LLT), named after Édouard Lucas and Derrick Henry Lehmer, is used for deterministically proving the primality of Mersenne numbers `M_p = 2^p - 1`, stating that:

`M_p` is prime if and only if `V_(M_p)(4,1) ≡ 4 mod M_p`

where `V_n(P,Q)` is the Lucas V sequence.

Taking advantage of the fact that `M_p + 1` is a power of `2`, and the Jacobi symbol `(12/(M_p)) = -1`, the test can be implemented more efficiently using the following iteration:

`S_0 = 4`; `S_n = (S_(n-1)^2 - 2) mod M_p`

then `M_p` is prime if and only if:

`S_(p-2) ≡ 0 mod M_p`

# Pocklington `n-1` primality proving

The Pocklington-Lehmer primality test, named after Henry Cabourn Pocklington and Derrick Henry Lehmer, is an improved version of the Lucas n-1 primality test, proving the primality of `n`, requiring only the partial factorization of `n-1`.

The basic version of this test states that, if `n > 1` is a positive odd integer, and there exists an integer `a` and a prime `p`, where `p | (n-1)` and `p > floor(sqrt(n))`, such that:

`a^(n-1) ≡ 1 mod n`

and

`gcd(a^((n-1)/p) - 1, n) = 1`

then `n` is prime.

In the generalized version, we write `n-1 = A * B`, where `gcd(A,B) = 1` and `A > B`, with the prime factorization of `A` known.

If for each prime factor `p` of `A` there exists an integer `a_p` such that:

`a_p^(n-1) ≡ 1 mod n`

and

`gcd(a_p^((n-1)/p)-1, n) = 1`

then `n` is prime.

The Fermat test can also be replaced with a Miller-Rabin test.

The beauty of this algorithm, is that we can apply it recursively on the prime factors of `n-1`, until we hit the base case of `2`. Doing this, we produce a primality certificate.

In practice, we can set the base case to `2^64`, where for inputs less than this bound, the Baillie-PSW primality test is known to be 100% correct.

In factoring `n-1` into `n-1 = A * B`, we can use trial division over some small primes, followed by the Lenstra elliptic-curve factorization method (ECM), finding one factor at a time, until we have `A > sqrt(n)` with all the prime factors of `A` known, recursively testing each prime factor of `A` for primality.

# Lucas `n+1` primality proving

Based on the ideas from the Pocklington `n-1` test, we can use the modular Lucas sequence `U_(n+1)(P,Q) mod n` with discriminant `D = P^2 - 4Q` such that the Jacobi symbol `(D/n) = -1`, allowing us to recursively prove the primality of `n` using the partial factorization of `n+1`.

Let `n + 1 = A * B`, where `A > B` and `gcd(A,B) = 1` and the prime factorization of `A` is known.

Choose `P` and `Q` such that `D = P^2 - 4Q` is not a square modulo `n`. If there exists a Lucas sequence `U_(n+1)` of discriminant `D`, such that:

`U_(n+1) ≡ 0 mod n`

and for every prime factor `p` of `A`:

`gcd(U_((n+1)/p) , n) = 1`

then `n` is prime.

If no such sequence exists for a chosen `P` and `Q`, a new `P'` and `Q'` with the same `D` can be computed as:

`P' = P+2`
`Q' = P + Q + 1`

Note: the same discriminant `D` must be used for all prime factors `p` of `A`.

# Lucas-Pocklington `n±1` primality proving

The Pocklington `n-1` and the Lucas `n+1` methods, can be combined into a single one, recursively factoring either `n-1` or `n+1`, whichever is easier to factorize first.

We start by writing:

`n - 1 = A_1 B_1`

`n + 1 = A_2 B_2`

for some `A_1 >= 2` and `A_2 >= 2` (using trial division to find `A_1` and `A_2`).

Next, we try to move more factors from `B_1` into `A_1` and/or from `B_2` into `A_2`, until we have either `A_1 > B_1` or `A_2 > B_2`.

We can do this by using the Lenstra elliptic-curve factorization method (ECM) on the product `B_1 B_2`, which finds a factor of either `n-1` or `n+1`.

If we find `A_1 > B_1`, then we run the Pocklington `n-1` test.
If we find `A_2 > B_2`, then we run the Lucas `n+1` test.

Same as in the Pocklington `n-1` test, we can use the Baillie-PSW test for inputs below `2^64` as a base case.