Lucas–Lehmer primality test

In mathematics, the Lucas–Lehmer test (LLT) is a primality test for Mersenne numbers. The test was originally developed by Édouard Lucas in 1878^[1] and subsequently proved by Derrick Henry Lehmer in 1930.

The Lucas–Lehmer test works as follows. Let M_p = 2^p − 1 be the Mersenne number to test with p an odd prime. The primality of p can be efficiently checked with a simple algorithm like trial division since p is exponentially smaller than M_p. Define a sequence ${\displaystyle \{s_{i}\}}$ for all i ≥ 0 by

${\displaystyle s_{i}={\begin{cases}4&{\text{if }}i=0;\\s_{i-1}^{2}-2&{\text{otherwise.}}\end{cases}}}$

The first few terms of this sequence are 4, 14, 194, 37634, ... (sequence A003010 in the OEIS). Then M_p is prime if and only if

${\displaystyle s_{p-2}\equiv 0{\pmod {M_{p}}}.}$

The number s_p − 2 mod M_p is called the Lucas–Lehmer residue of p. (Some authors equivalently set s₁ = 4 and test s_p−1 mod M_p). In pseudocode, the test might be written as

// Determine if Mp = 2p − 1 is prime for p > 2
Lucas–Lehmer(p)
    var s = 4
    var M = 2p − 1
    repeat p − 2 times:
        s = ((s × s) − 2) mod M
    if s == 0 return PRIME else return COMPOSITE

Performing the mod M at each iteration ensures that all intermediate results are at most p bits (otherwise the number of bits would double each iteration). The same strategy is used in modular exponentiation.

Starting values s₀ other than 4 are possible, for instance 10, 52, and others (sequence A018844 in the OEIS).^[2] The Lucas-Lehmer residue calculated with these alternative starting values will still be zero if M_p is a Mersenne prime. However, the terms of the sequence will be different and a non-zero Lucas-Lehmer residue for non-prime M_p will have a different numerical value from the non-zero value calculated when s₀ = 4.

It is also possible to use the starting value (2 mod M_p)(3 mod M_p)⁻¹, usually denoted by 2/3 for short.^[2] This starting value equals (2^p + 1) /3, the Wagstaff number with exponent p.

Starting values like 4, 10, and 2/3 are universal, that is, they are valid for all (or nearly all) p. There are infinitely many additional universal starting values.^[2] However, some other starting values are only valid for a subset of all possible p, for example s₀ = 3 can be used if p = 3 (mod 4).^[3] This starting value was often used where suitable in the era of hand computation, including by Lucas in proving M₁₂₇ prime.^[4] The first few terms of the sequence are 3, 7, 47, ... (sequence A001566 in the OEIS).

If s_p−2 = 0 mod M_p then the penultimate term is s_p−3 = ± 2^(p+1)/2 mod M_p. The sign of this penultimate term is called the Lehmer symbol ϵ(s₀, p).

In 2000 S.Y. Gebre-Egziabher proved that for the starting value 2/3 and for p ≠ 5 the sign is:

${\displaystyle \epsilon ({2 \over 3},\ p)=(-1)^{p-1 \over 2}}$

That is, ϵ(2/3, p) = +1 if p = 1 (mod 4) and p ≠ 5.^[2]

The same author also proved Woltman's conjecture^[5] that the Lehmer symbols for starting values 4 and 10 when p is not 2 or 5 are related by:

${\displaystyle \epsilon (10,\ p)=\epsilon (4,\ p)\ \times \ (-1)^{{(p+1)(p+3)} \over 8}}$

That is, ϵ(4, p) × ϵ(10, p) = 1 if p = 5 or 7 (mod 8) and p ≠ 2, 5.^[2]

OEIS sequence A123271 shows ϵ(4, p) for each Mersenne prime M_p.

In the algorithm as written above, there are two expensive operations during each iteration: the multiplication s × s, and the mod M operation. The mod M operation can be made particularly efficient on standard binary computers by observing that

${\displaystyle k\equiv (k\,{\bmod {\,}}2^{n})+\lfloor k/2^{n}\rfloor {\pmod {2^{n}-1}}.}$

This says that the least significant n bits of k plus the remaining bits of k are equivalent to k modulo 2ⁿ−1. This equivalence can be used repeatedly until at most n bits remain. In this way, the remainder after dividing k by the Mersenne number 2ⁿ−1 is computed without using division. For example,

Moreover, since s × s will never exceed M² < 2^2p, this simple technique converges in at most 1 p-bit addition (and possibly a carry from the pth bit to the 1st bit), which can be done in linear time. This algorithm has a small exceptional case. It will produce 2ⁿ−1 for a multiple of the modulus rather than the correct value of 0. However, this case is easy to detect and correct.

With the modulus out of the way, the asymptotic complexity of the algorithm only depends on the multiplication algorithm used to square s at each step. The simple "grade-school" algorithm for multiplication requires O(p²) bit-level or word-level operations to square a p-bit number. Since this happens O(p) times, the total time complexity is O(p³). A more efficient multiplication algorithm is the Schönhage–Strassen algorithm, which is based on the Fast Fourier transform. It only requires O(p log p log log p) time to square a p-bit number. This reduces the complexity to O(p² log p log log p) or Õ(p²).^[6] An even more efficient multiplication algorithm, Fürer's algorithm, only needs ${\displaystyle p\log p\ 2^{O(\log ^{*}p)}}$ time to multiply two p-bit numbers.

By comparison, the most efficient randomized primality test for general integers, the Miller–Rabin primality test, requires O(k n² log n log log n)^{[contradictory]} bit operations using FFT multiplication for an n-digit number, where k is the number of iterations and is related to the error rate. For constant k, this is in the same complexity class as the Lucas-Lehmer test. In practice however, the cost of doing many iterations and other differences leads to worse performance for Miller–Rabin.^{[clarification needed]} The most efficient deterministic primality test for any n-digit number, the AKS primality test, requires Õ(n⁶) bit operations in its best known variant and is extremely slow even for relatively small values.

The Mersenne number M₃ = 2³−1 = 7 is prime. The Lucas–Lehmer test verifies this as follows. Initially s is set to 4 and then is updated 3−2 = 1 time:

• s ← ((4 × 4) − 2) mod 7 = 0.

Since the final value of s is 0, the conclusion is that M₃ is prime.

On the other hand, M₁₁ = 2047 = 23 × 89 is not prime. Again, s is set to 4 but is now updated 11−2 = 9 times:

• s ← ((4 × 4) − 2) mod 2047 = 14
• s ← ((14 × 14) − 2) mod 2047 = 194
• s ← ((194 × 194) − 2) mod 2047 = 788
• s ← ((788 × 788) − 2) mod 2047 = 701
• s ← ((701 × 701) − 2) mod 2047 = 119
• s ← ((119 × 119) − 2) mod 2047 = 1877
• s ← ((1877 × 1877) − 2) mod 2047 = 240
• s ← ((240 × 240) − 2) mod 2047 = 282
• s ← ((282 × 282) − 2) mod 2047 = 1736

Since the final value of s is not 0, the conclusion is that M₁₁ = 2047 is not prime. Although M₁₁ = 2047 has nontrivial factors, the Lucas–Lehmer test gives no indication about what they might be.

The proof of correctness for this test presented here is simpler than the original proof given by Lehmer. Recall the definition

${\displaystyle s_{i}={\begin{cases}4&{\text{if }}i=0;\\s_{i-1}^{2}-2&{\text{otherwise.}}\end{cases}}}$

The goal is to show that M_p is prime iff ${\displaystyle s_{p-2}\equiv 0{\pmod {M_{p}}}.}$

The sequence ${\displaystyle {\langle }s_{i}{\rangle }}$ is a recurrence relation with a closed-form solution. Let ${\displaystyle \omega =2+{\sqrt {3}}}$ and ${\displaystyle {\bar {\omega }}=2-{\sqrt {3}}}$. It then follows by induction that ${\displaystyle s_{i}=\omega ^{2^{i}}+{\bar {\omega }}^{2^{i}}}$ for all i:

${\displaystyle s_{0}=\omega ^{2^{0}}+{\bar {\omega }}^{2^{0}}=\left(2+{\sqrt {3}}\right)+\left(2-{\sqrt {3}}\right)=4}$

and

${\displaystyle {\begin{aligned}s_{n}&=s_{n-1}^{2}-2\\&=\left(\omega ^{2^{n-1}}+{\bar {\omega }}^{2^{n-1}}\right)^{2}-2\\&=\omega ^{2^{n}}+{\bar {\omega }}^{2^{n}}+2(\omega {\bar {\omega }})^{2^{n-1}}-2\\&=\omega ^{2^{n}}+{\bar {\omega }}^{2^{n}}.\end{aligned}}}$

The last step uses ${\displaystyle \omega {\bar {\omega }}=\left(2+{\sqrt {3}}\right)\left(2-{\sqrt {3}}\right)=1.}$ This closed form is used in both the proof of sufficiency and the proof of necessity.

The goal is to show that ${\displaystyle s_{p-2}\equiv 0{\pmod {M_{p}}}}$ implies that ${\displaystyle M_{p}}$ is prime. What follows is a straightforward proof exploiting elementary group theory given by J. W. Bruce^[7] as related by Jason Wojciechowski.^[8]

Suppose ${\displaystyle s_{p-2}\equiv 0{\pmod {M_{p}}}.}$ Then

${\displaystyle \omega ^{2^{p-2}}+{\bar {\omega }}^{2^{p-2}}=kM_{p}}$

for some integer k, so

${\displaystyle \omega ^{2^{p-2}}=kM_{p}-{\bar {\omega }}^{2^{p-2}}.}$

Multiplying by ${\displaystyle \omega ^{2^{p-2}}}$ gives

${\displaystyle \left(\omega ^{2^{p-2}}\right)^{2}=kM_{p}\omega ^{2^{p-2}}-(\omega {\bar {\omega }})^{2^{p-2}}.}$

Thus,

${\displaystyle \omega ^{2^{p-1}}=kM_{p}\omega ^{2^{p-2}}-1.\qquad \qquad (1)}$

For a contradiction, suppose M_p is composite, and let q be the smallest prime factor of M_p. Mersenne numbers are odd, so q > 2. Let ${\displaystyle \mathbb {Z} _{q}}$ be the integers modulo q, and let ${\displaystyle X=\left\{a+b{\sqrt {3}}\mid a,b\in \mathbb {Z} _{q}\right\}.}$ Multiplication in ${\displaystyle X}$ is defined as

${\displaystyle \left(a+{\sqrt {3}}b\right)\left(c+{\sqrt {3}}d\right)=[(ac+3bd)\,{\bmod {\,}}q]+{\sqrt {3}}[(ad+bc)\,{\bmod {\,}}q].}$^{[note 1]}

Clearly, this multiplication is closed, i.e. the product of numbers from X is itself in X. The size of X is denoted by ${\displaystyle |X|.}$

Since q > 2, it follows that ${\displaystyle \omega }$ and ${\displaystyle {\bar {\omega }}}$ are in X.^{[note 2]} The subset of elements in X with inverses forms a group, which is denoted by X* and has size ${\displaystyle |X^{*}|.}$ One element in X that does not have an inverse is 0, so ${\displaystyle |X^{*}|\leq |X|-1=q^{2}-1.}$

Now ${\displaystyle M_{p}\equiv 0{\pmod {q}}}$ and ${\displaystyle \omega \in X}$, so

${\displaystyle kM_{p}\omega ^{2^{p-2}}=0}$

in X. Then by equation (1),

${\displaystyle \omega ^{2^{p-1}}=-1}$

in X, and squaring both sides gives

${\displaystyle \omega ^{2^{p}}=1.}$

Thus ${\displaystyle \omega }$ lies in X* and has inverse ${\displaystyle \omega ^{2^{p}-1}.}$ Furthermore, the order of ${\displaystyle \omega }$ divides ${\displaystyle 2^{p}.}$ However ${\displaystyle \omega ^{2^{p-1}}\neq 1}$, so the order does not divide ${\displaystyle 2^{p-1}.}$ Thus, the order is exactly ${\displaystyle 2^{p}.}$

The order of an element is at most the order (size) of the group, so

${\displaystyle 2^{p}\leq |X^{*}|\leq q^{2}-1<q^{2}.}$

But q is the smallest prime factor of the composite ${\displaystyle M_{p}}$, so

${\displaystyle q^{2}\leq M_{p}=2^{p}-1.}$

This yields the contradiction ${\displaystyle 2^{p}<2^{p}-1}$. Therefore, ${\displaystyle M_{p}}$ is prime.

In the other direction, the goal is to show that the primality of ${\displaystyle M_{p}}$ implies ${\displaystyle s_{p-2}\equiv 0{\pmod {M_{p}}}}$. The following simplified proof is by Öystein J. Rödseth.^[9]

Since ${\displaystyle 2^{p}-1\equiv 7{\pmod {12}}}$ for odd ${\displaystyle p>1}$, it follows from properties of the Legendre symbol that ${\displaystyle (3|M_{p})=-1.}$ This means that 3 is a quadratic nonresidue modulo ${\displaystyle M_{p}.}$ By Euler's criterion, this is equivalent to

${\displaystyle 3^{\frac {M_{p}-1}{2}}\equiv -1{\pmod {M_{p}}}.}$

In contrast, 2 is a quadratic residue modulo ${\displaystyle M_{p}}$ since ${\displaystyle 2^{p}\equiv 1{\pmod {M_{p}}}}$ and so ${\displaystyle 2\equiv 2^{p+1}=\left(2^{\frac {p+1}{2}}\right)^{2}{\pmod {M_{p}}}.}$ This time, Euler's criterion gives

${\displaystyle 2^{\frac {M_{p}-1}{2}}\equiv 1{\pmod {M_{p}}}.}$

Combining these two equivalence relations yields

${\displaystyle 24^{\frac {M_{p}-1}{2}}\equiv \left(2^{\frac {M_{p}-1}{2}}\right)^{3}\left(3^{\frac {M_{p}-1}{2}}\right)\equiv (1)^{3}(-1)\equiv -1{\pmod {M_{p}}}.}$

Let ${\displaystyle \sigma =2{\sqrt {3}}}$, and define X as before as the ring ${\displaystyle X=\{a+b{\sqrt {3}}\mid a,b\in \mathbb {Z} _{M_{p}}\}.}$ Then in the ring X, it follows that

${\displaystyle {\begin{aligned}(6+\sigma )^{M_{p}}&=6^{M_{p}}+\left(2^{M_{p}}\right)\left({\sqrt {3}}^{M_{p}}\right)\\&=6+2\left(3^{\frac {M_{p}-1}{2}}\right){\sqrt {3}}\\&=6+2(-1){\sqrt {3}}\\&=6-\sigma ,\end{aligned}}}$

where the first equality uses the Binomial Theorem in a finite field, which is

${\displaystyle (x+y)^{M_{p}}\equiv x^{M_{p}}+y^{M_{p}}{\pmod {M_{p}}}}$,

and the second equality uses Fermat's little theorem, which is

${\displaystyle a^{M_{p}}\equiv a{\pmod {M_{p}}}}$

for any integer a. The value of ${\displaystyle \sigma }$ was chosen so that ${\displaystyle \omega ={\frac {(6+\sigma )^{2}}{24}}.}$ Consequently, this can be used to compute ${\displaystyle \omega ^{\frac {M_{p}+1}{2}}}$ in the ring X as

${\displaystyle {\begin{aligned}\omega ^{\frac {M_{p}+1}{2}}&={\frac {(6+\sigma )^{M_{p}+1}}{24^{\frac {M_{p}+1}{2}}}}\\&={\frac {(6+\sigma )(6+\sigma )^{M_{p}}}{24\cdot 24^{\frac {M_{p}-1}{2}}}}\\&={\frac {(6+\sigma )(6-\sigma )}{-24}}\\&=-1.\end{aligned}}}$

All that remains is to multiply both sides of this equation by ${\displaystyle {\bar {\omega }}^{\frac {M_{p}+1}{4}}}$ and use ${\displaystyle \omega {\bar {\omega }}=1}$, which gives

${\displaystyle {\begin{aligned}\omega ^{\frac {M_{p}+1}{2}}{\bar {\omega }}^{\frac {M_{p}+1}{4}}&=-{\bar {\omega }}^{\frac {M_{p}+1}{4}}\\\omega ^{\frac {M_{p}+1}{4}}+{\bar {\omega }}^{\frac {M_{p}+1}{4}}&=0\\\omega ^{\frac {2^{p}-1+1}{4}}+{\bar {\omega }}^{\frac {2^{p}-1+1}{4}}&=0\\\omega ^{2^{p-2}}+{\bar {\omega }}^{2^{p-2}}&=0\\s_{p-2}&=0.\end{aligned}}}$

Since ${\displaystyle s_{p-2}}$ is 0 in X, it is also 0 modulo ${\displaystyle M_{p}.}$

The Lucas–Lehmer test is one of the main primality tests used by the Great Internet Mersenne Prime Search (GIMPS) to locate large primes. This search has been successful in locating many of the largest primes known to date.^[10] The test is considered valuable because it can provably test a large set of very large numbers for primality within an affordable amount of time. In contrast, the equivalently fast Pépin's test for any Fermat number can only be used on a much smaller set of very large numbers before reaching computational limits.

• Mersenne's conjecture
• Lucas–Lehmer–Riesel test

• Crandall, Richard; Pomerance, Carl (2001), "Section 4.2.1: The Lucas–Lehmer test", Prime Numbers: A Computational Perspective (1st ed.), Berlin: Springer, pp. 167–170, ISBN 0-387-94777-9

• Weisstein, Eric W. "Lucas–Lehmer test". MathWorld.
• GIMPS (The Great Internet Mersenne Prime Search)
• A proof of Lucas–Lehmer–Reix test (for Fermat numbers)
• Lucas–Lehmer test at MersenneWiki