### 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:

1. a^d ≡ 1 mod n
2. 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:
1. Perform a Miller-Rabin test to base 2.
2. Find the first P >= 1 for which the Jacobi symbol ((P^2+4)/n) = -1.
3. 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 PomeranceJohn Selfridge, and Samuel Wagstaff, is a combination of two primality tests:
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):
1. U_d(P,Q) ≡ 0 mod n
2. 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):
1. U_d(P,Q) ≡ 0 mod n and V_d(P,Q) ≡ ±2 mod n
2. 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:
1. find an integer a such that the Jacobi symbol (a/n) = -1.
2. 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.

### # Implementations

This section includes source-code implementations of the algorithms described in this post.

:: Sidef
:: Perl

Additioanlly, some bits and pieces used in primality testing algorithms:

:: Sidef
:: Perl