理论

设 $D$ 为无平方整数，记 $O$ 为二次域 $\mathbb{Q}(\sqrt{D})$ 上的 代数整数环，有类比的定理：

【二次域中的 Fermat 小定理】设 $a \in O$ ， $p$ 为奇素数，且 $p \nmid d$ ，则

$$ a^p = \begin{cases} a, &\text{若 } \left(\frac{D}{p}\right) = 1 \\ \overline{a}, &\text{若 } \left(\frac{D}{p}\right) = -1\end{cases} \pmod p $$

由于乘幂并不容易计算，引入一类特殊的二阶递推序列：设 $P, Q \in \mathbb{Z}$ ，特征方程 $x^2 - Px+ Q = 0$ 的两根分别为 $a, b$ ，定义 $P, Q$ 对应的 Lucas 序列为

$$ U_n = \frac{a^n - b^n}{a - b},\quad V_n = a^n + b^n $$

我没理解乘幂计算的在哪里有困难，不能快速幂吗？猜测是那个除 $2$ 比较头疼。

其与二次域的关系是

$$ \left(\frac{P + \sqrt{D}}{2}\right)^n = \frac{V_n + U_n \sqrt{D}}{2}, \quad D = P^2 - 4Q $$

首先我们限定 $\gcd(QD, p) = 1$ ，保证 $a,b,a-b$ 在模 $p$ 下可逆。再记 $\varepsilon = \left(\frac{D}{p}\right)$ ，不难根据 Fermat 小定理证明 $U_{p - \varepsilon} \equiv 0 \pmod p$ 。

算法描述

输入奇数 $N$ ，选取 $P, Q$ 作为 Lucas 序列参数，使得 $\varepsilon = \left(\frac{D}{N}\right) \neq 0$ 。

1. 若 $U_{N - \varepsilon} \equiv 0 \pmod N$ ，则输出 $N$ 可能为奇素数。
2. 若 $U_{N - \varepsilon} \not\equiv 0 \pmod N$ ，则输出 $N$ 为合数。

由于 $2$ 在奇数下的逆元总是存在的，我们可以使用矩阵快速幂计算上式。

对于 $P,Q$ 的选取，有两种建议：

1. Selfridge 建议：取 $D$ 为 $5, -7, 9, -11,13,\cdots$ 中使得 $\left(\frac{D}{N}\right) = -1$ 的第一个数，选择 $P = 1, Q = \frac{1-D}{2}$ 。
2. Baillie 建议：取 $D$ 为 $5, 9, 13, 17, 21, \cdots$ 中使得使得 $\left(\frac{D}{N}\right) = -1$ 的第一个数，选择 $P$ 为大于 $\sqrt{D}$ 的最小奇数， $Q = \frac{P^2 - D}{4}$ 。

论文中提到，取到合适的 $D$ 的期望次数为 $2$ 。若尝试次数过多，则说明 $N$ 很可能是平方数，需要检测下。

备注

网上找到的很多 Lucas 检测代码都是关于 $n+1$ 的因式分解，我要是能分解还判素数干嘛。。

似乎大家都建议固定 $D$ ，然后 $(p, q) \to (p + 2, p + q + 1)$ 这样，但是通不过 LOJ143，而且反例很多的样子。。。反正计算 Jacobi 的复杂度是 $O(\log^2 n)$ 的，不如换 $D$ 。

参考

• 素数判定 - mathmu
• 如何编程判断一个数是否是质数？ - Wiener ES的回答
• Primality Proving 3.2 n+1 tests and the Lucas-Lehmer test

代码

#ifndef ALGO_MATH_LUCAS_TEST
#define ALGO_MATH_LUCAS_TEST

#include "../../base.hpp"
#include "../qpow/u64.hpp"
#include "../jacobi.hpp"
#include "../square-test.hpp"

#include <functional>

bool lucas_test_pd(u64 N, u64 P, u64 D) { // jacobi(D, N) == -1
  using pii = std::pair<u64, u64>;
  auto shift2 = [N](u64 n) {
    return n % 2 == 0 ? n / 2 : (N + n) / 2;
  };
  std::function<pii(u64)> calc = [&](u64 n) {
    if (n == 1)
      return pii{P, 1};
    auto [v, u] = calc(n / 2);
    u64 u2 = u128(u) * v % N;
    u64 v2 = (u128(v) * v + u128(u) * u % N * D) % N;
    v2 = shift2(v2);
    if (n % 2 == 1) {
      u64 u21 = (u128(P) * u2 + v2) % N;
      u64 v21 = (u128(D) * u2 + u128(P) * v2) % N;
      u2 = shift2(u21), v2 = shift2(v21);
    }
    return pii{v2, u2};
  };
  auto [v, u] = calc(N + 1);
  return u == 0;
}

bool lucas_test(u64 n) {
  if (n <= 6)
    return n == 2 || n == 3 || n == 5;
  if (n % 6 != 1 && n % 6 != 5)
    return false;
  if (square_test(n))
    return false;
  for (u64 d = 5, c = 0; d < n && c < 4; d += 4) {
    i32 j = jacobi(d, n);
    if (j == 0)
      return false;
    if (j == 1)
      continue;
    if (!lucas_test_pd(n, d, d))
      return false;
    c++;
  }
  return true;
}

#endif // ALGO_MATH_LUCAS_TEST