有报告称，新的梅森素数即将被发现。（已经发现了，是 2^{136279841}-1）
那么梅森数如何进行素性检验呢？AKS 算法虽然快，但是对于超大的梅森数也显得捉襟见肘。对此我们有更好的检验方法：Lucas–Lehmer 检验。
设数列 s_0=4，s_i=s_{i-1}^2-2。则对于奇素数 p，M_p=2^p-1 是素数当且仅当 M_p\mid s_{p-2}。
这个方法惊人地简单，而且其背后的证明也并不复杂，我们接下来做详细介绍。
在此之前，我们尝试得出 s_n 的封闭形式，设 \omega=2+\sqrt 3，\bar \omega=2-\sqrt 3，则 s_n=\omega^{2^n}+{\bar \omega}^{2^n}。这可以简单地归纳证明（首先 s_0 符合归纳假设）：
\begin{aligned}s_{n+1}&=s_n^2-2 \\ &=\omega^{2^{n+1}}+{\bar \omega}^{2^{n+1}}+2(\omega \bar \omega)^{2^n}-2 \\ &=\omega^{2^{n+1}}+{\bar \omega}^{2^{n+1}} \end{aligned}
最后一步是因为 \omega \bar\omega=1。这个封闭形式在后面的证明中很重要。
充分性的证明
现在我们有 M_p\mid s_{p-2}，可以设 s_{p-2}=kM_p，即
\omega^{2^{p-2}}+{\bar \omega}^{2^{p-2}}=kM_p
等式两侧同乘 \omega^{2^{p-2}} 就有
\omega^{2^{p-1}}+1=kM_p \omega^{2^{p-2}}
考虑反证法，假设 M_p 为合数，并设 q 为其最小质因子（显然会有 q>2）。然后我们考虑这样的一个代数整数环 R：其中所有元素都是 a+b\sqrt 3 的形式（0\leq a,b < q 且 a,b 为整数），这个环上的乘法定义为
(a+b\sqrt 3)(c+d\sqrt 3)=((ac+3bd)\bmod q)+((ad+bc)\bmod q)\sqrt 3
由于 q\mid M_p，所以在 R 上有
\omega^{2^{p-1}}=-1
等式两边同时平方，就有 \omega^{2^p}=1，这意味着 \omega 的阶就是 2^p，不可能更小（由于任意 t\mid 2^p 都有 \omega^t \neq 1，若不然会与 \omega^{2^{p-1}}=-1 相悖）。
另一方面，R 中任意元素的阶都不可能超过 q^2，即 R 中元素的数量 （对于 x\in R，连接一条有向边 x\to x\omega，从 1 开始最多走 q^2 步回到 1）。
于是就得到了关系：
\text{ord}(\omega)=2^p\leq q^2\leq M_p=2^p-1
这样就导出了矛盾，故 M_p 必须为素数。
必要性的证明
现在我们已知的条件是 M_p 为素数。
对于奇数 p>1 有 2^p-1\equiv 7 \pmod{12}，根据二次互反律（下面用到的 Legendre 符号在此文开头有介绍）：
\left( \frac{3}{M_p}\right)\left( \frac{M_p}{3}\right) =(-1)^{\frac{M_p-1}{2}\frac{3-1}{2}}=-1
而
\left( \frac{M_p}{3}\right)\equiv M_p\pmod 3
显然 M_p\bmod 3 = 1，所以 \left( \dfrac{3}{M_p}\right)=-1，即
3^{\frac{M_p-1}{2}}\equiv -1 \pmod{M_p}
由于 2^p\equiv 1 \pmod{M_p}，我们有
2^{\frac{M_p-1}{2}}\equiv 2^{(2^{p-1}-1)\bmod p}\pmod{M_p}
其中 2^{p-1}\bmod p=1，这是费马小定理，于是
2^{\frac{M_p-1}{2}}\equiv 1 \pmod{M_p}
根据这两条结果，我们先来尝试计算 \omega^{2^{p-1}}\bmod M_p。为了方便表示，类比构建充分性证明中的环 R，只不过模数从 q 改为了 M_p。自然下面的运算都在 R 上考虑。
设 \sigma = 2\sqrt 3，则 \omega=\dfrac{(6+\sigma)^2}{24}：
\omega^{2^{p-1}}=\omega^{\frac{M_p+1}{2}}=\frac{(6+\sigma)^{M_p+1}}{24\times 24^{\frac{M_p-1}{2}}}
分子和分母可以分开来看，分母上有 24=2^3\times 3，故 24^{\frac{M_p-1}{2}}= -1。
分子上有
\begin{aligned}(6+\sigma)^{M_p}&= 6^{M_p}+2^{M_p}3^{\frac{M_p-1}{2}}\sqrt 3 \\ &= 6-2\sqrt 3 \end{aligned}
所以
\omega^{2^{p-1}}= \frac{(6+\sigma)(6-\sigma)}{-24}=-1
对等式两侧同乘以 {\bar \omega}^{2^{p-2}} 并利用 \omega {\bar \omega}=1，就得到了
\omega^{2^{p-2}}+{\bar \omega}^{2^{p-2}}=0
这也就证明了必要性。
一些补充内容
首先是具体实现上的优化，考虑如何快速计算 n \bmod M_p。根据
n=(n \bmod 2^p)+2^p\lfloor n/2^p\rfloor
等式两边同时对 M_p 取模，就得到了
n \equiv (n\bmod 2^p)+\lfloor n/2^p \rfloor \pmod{M_p}
这在 n 用二进制表示时是很容易操作的。
而且由于计算中每一项都要取模，所以实操中 n <M_p^2，只用 \Theta(p) 的时间就可以完成取模。
总计算时间复杂度还是 \Theta(p^2 \log p) 的。
其次，真的没办法再优化了吗？
如果还能优化，那么人类就又在梅森素数的探索上前进了一大步。如果是对于一般的递推，那就更为困难了。要知道著名的 Mandelbrot 集合就是根据递推 a_n=a_{n-1}^2+C 的模是否发散来生成的（若不发散到无穷，则 C 属于该集合）。