### （二）、Lucas 序列素性测试
若 $P \in \mathbb{Z}+，Q \in \mathbb{Z}，D=P^2-4Q$为整数，定义 Lucas 序列 $U，V$：
$U_0=0，U_1=1，\dots，U_{k+2}=PU_{k+1}-QU_k （k \ge 0）$，
$V_0=2，V_1=P，\dots，V_{k+2}=PV_{k+1}-QV_k （k \ge 0）$，
若 $n \in N+$ 是大于 1 的奇数，选择 $D，P，Q$ 满足 Jacobi 符号 $J(\frac{D}{n}) = -1$，
当 $n$ 是素数且 $GCD(n，Q)=1$，我们有：
$$U_{n+1} = 0 (mod\ \ n)\tag{5-1}$$ $$V_{n+1} = 2Q (mod\ \ n)\tag{5-2}$$
另外，如果 $n$ 是奇素数，$\delta(n) = n - J(\frac{D}{n})$，且 $GCD(D, n) = 1$，则：
$$U_{\delta(n)} \equiv 0 (mod\ \ n)\tag{5-3}$$ $$V_{\delta(n)} \equiv 2Q^{(1-J(\frac{D}{n}))/2} (mod\ \ n)\tag{5-4}$$ $$U_n \equiv J(\frac{D}{n}) (mod\ \ n)\tag{5-5}$$ $$V_n \equiv V_1 = P (mod\ \ n)\tag{5-6}$$
利用上面的结果的素性测试，称为 Lucas 序列素性测试。如果合数通过测试 $(5-1)$，称为以 $P，Q$ 为参数的 lucas 伪素数，简称 lpsp(P，Q)。如果合数通过测试 $(5-2)$，称为以 $P，Q$ 为参数的 lucas-v 伪素数，简称 vpsp(P，Q)。
对 Lucas 序列有如下快速迭代公式。
$$U_{2k} = U_kV_k\tag{5-7}$$ $$V_{2k} = V_k^2 - 2Q^k\tag{5-8}$$ $$U_{k+1} = (PU_k + V_k) / 2\tag{5-9}$$ $$V_{k+1} = (DU_k + PV_k) / 2\tag{5-10}$$
#### A、参数的选择
方法一、$D$ 在序列 $5， -7， 9， -11， 13， -15，\dots$ 中选择第一个满足 $J(\frac{D}{n}) = -1$ 的值，令 $P = 1, Q = \frac{1 - D}{4}$
这个方法，不会设置 $Q = 1$，但是常会设置 $Q = -1$，当 $n \equiv \pm 3 (mod\ \ 10)$ 时，此时 $D = 5$。
计算表明， 当 $Q \equiv \pm 1 (mod\ \ n)$ 时要比 $Q \not\equiv \pm 1 (mod\ \ n)$ 有更多合数满足公式 $(5-3)$、$(5-4)$、$(5-5)$、$(5-6)$，因此有下面的改进方法。
方法二、按照方法一的选择参数，当出现 $Q = -1$时候，置 $D，P，Q$ 都为5。
用方法二的参数，前 10 个 lpsp 是 323，377，1159，1829，3827，5459，5777，9071，9179，10877，$10^{15}$ 以内有 2402549 个 lpsp，vpsp 更稀少，$10^{15}$ 以内只有 5 个 vpsp。
#### B、强 lucas 素性测试
如果 $n$ 是奇素数，$J(\frac{D}{n}) = -1$，$n + 1 = q * 2^s$，$q$ 是奇数，则必有 $r$ 满足 $0 \leq r < s$ 满足下面两个条件之一： $$U_q \equiv 0 (mod\ \ n) \tag{5-11}$$ $$V_{q*2^r} \equiv 0 (mod\ \ n) \tag{5-12}$$
如果 $n$ 是合数，用方法二选择参数$D、P、Q$，$J(\frac{D}{n}) = -1$，满足 $(5-11)$ 或者 $(5-12)$，$n$ 称为强 lucas 伪素数，记做 slpsp(P，Q)。前 10 个 slpsp 是5459，5777，10877，16109，18971，22499，24569，25199，40309，58519。强 lucas 伪素数远少于 lucas 伪素数，$10^{15}$ 以内只有 474971 个。
#### C、BPSW 算法
该算法是对 $n$ 首先进行基 2 的 Miller-Rabin 测试，然后用方法二选择参数进行强 Lucas 素性测试。即：
1、如果 $n$ 没有通过基 2 Miller-Rabin 测试，输出 $n$ 是合数，结束。
2、用方法二选择参数 $D，P，Q$。
3、进行强 lucas 素性测试，通过则输出 $n$ 可能是素数，否则 $n$ 是合数。
#### D、增强的 BPSW 算法
基于 vpsp 很稀少的事实，针对 BPSW 增加了两个检验。
1、如果 $n$ 没有通过基 2 Miller-Rabin 测试，输出 $n$ 是合数，结束。
2、用方法二选择参数 $D，P，Q$。
3、进行强 lucas 素性测试，通过则转步骤 4 进一步检验，否则输出 $n$ 是合数，结束。
4、基于步骤 3 的计算，进一步计算并检验 $5-2$，不满足 $5-2$，则输出 $n$ 是合数，结束。
5、验证 $n$ 是不是满足 $Q^{(n+1)/2} \equiv Q*J(\frac{Q}{n}) (mod\ \ n)$，如果满足，输出 $n$ 可能素数，否则输出 $n$ 是合数。
#### E、实现算法：
输入待测试正奇数 n
1、d = 5, deta = 2
2、测试 $J(\frac{d}{n})$，如果不等于-1，则 deta = -deta, d = deta - d，重复 2。
3、置 D = d, P = 1, Q = (1 - D) / 4，如果 Q = -1，则 D = 5，P = 5，Q = 5。
4、令 $n + 1 = q * 2^s$，$q$ 是奇数，令 $qs$ 表示 $q$ 的二进制表示的长度，即 $2^{qs-1} \le q < 2^{qs}$。$U_c = U_1 = 1, V_c = V_1 = P, Q_c = Q$。
5、若$qs$ > 1，令 $i$ 从 $qs - 2$ 到 $0$ 循环执行 5.1-5.2（$qs$ 至少为1，且最高位肯定是1，所以这里从第二个二进制位开始）
5.1 $U_c = U_c * V_c (mod\ \ n)，V_c = V_c^2 - 2Q_c (mod\ \ n)，Q_c = Q_c^2 (mod\ \ n)$
5.2 如果 $q$ 二进制表示第 $i$ 位是 $1$，则 $U_c = (PU_c + V_c)/2 (mod\ \ n)，V_c = (DU_c + PV_c)/2 (mod\ \ n)，Q_c = Q_c Q (mod\ \ n) $。
6、此时 $U_c$ 即为 $U_q$，若 $U_c \not = {0}$，执行7，否则跳到9。
7、继续测试 $V_c$。令 $r = 0$，循环执行：
7.1、如果 $V_c = 0$，跳转到9。
7.2、$V_c = V_c^2 - 2Q_c (mod\ \ n)，Q_c = Q_c^2 (mod\ \ n)， r = r + 1$。如果 $r = s$，跳转到8，否则跳转到 7.1 继续循环。
8、输出 $n$ 是合数，结束。
9、此时，强 lucas 素性测试输出 $n$ 可能是素数，如果进一步测试，进行步骤 10。
10、强化的 BPSW 算法，利用了 5-7 的步骤的计算，附加计算并不多。分成两个检验：
10.1、循环执行 $V_c = V_c^2 - 2Q_c (mod\ \ n)，Q_c = Q_c^2 (mod\ \ n)， r = r + 1$，直到 $r = s$。
10.2、此时，$V_c = V_{n+1}$，测试 $V_c \equiv 2Q (mod\ \ n)$ 是否成立，如果不成立，输出 $n$ 是合数，结束。
10.3、测试 $Q^{(n+1)/2} \equiv Q*J(\frac{Q}{n}) (mod\ \ n)$ 是否成立，如果不成立，输出 $n$ 是合数，结束。否则，输出 $n$ 通过强化的 BPSW 素性测试，可能是素数。
mathematica 发表于 2019-9-29 09:19
我觉得你可真够无聊的,其实
baillie psw算法已经够好的了,
https://bbs.emath.ac.cn/forum.php?mod=vie ...
刚写完BPSW算法，发现你吹的BPSW照样对 1093^2 和 3511^2 无能，不知道为什么文献故意忽略掉这事~ @无心人 我在前面已经把完全平方数排除了，怎么无能了？？？？？？？？？
nyy 发表于 2022-8-30 12:44
刚写完BPSW算法，发现你吹的BPSW照样对 1093^2 和 3511^2 无能，不知道为什么文献故意忽略掉这事~ @无心 ...
标准BPSW只做基2米勒罗宾测试，这俩平方和含这俩平方因子的的一些数，滤不掉，后续的，求Jacobi符号的算法，对有平方因子的会持续失败
无心人 发表于 2022-8-31 16:03
标准BPSW只做基2米勒罗宾测试，这俩平方和含这俩平方因子的的一些数，滤不掉，后续的，求Jacobi符号的算 ...
Rabin-Miller 强伪素数可能含平方因子吗？
https://bbs.emath.ac.cn/forum.php?mod=viewthread&tid=3710&fromuid=14149
(出处: 数学研发论坛)
你和郭先强一样，都没搞懂！
看看这个数
26092328809=1093*1093*21841
{{1093, 2}, {21841, 1}}
计算雅克比符号
JacobiSymbol
得到-1
质因数分解
5654273717={{1093, 2}, {4733, 1}}
计算雅克比符号
JacobiSymbol
-1
最保守的办法，就是开平方，上面两个例子证明：含有平方因子，根本就不影响。
还有就是你可以先计算一万次雅克比符号，如果没获得雅克比符号等于-1的，那么赶快就开一次平方看看，
或者你一开始就开平方看看。
或者你直接除以1093、3511这两个素数，
再看601401837037=1093*1093*503413
质因数分解，得到{{1093, 2}, {503413, 1}}
计算雅克比符号：JacobiSymbol
得到-1，
这都证明是否含有完全平方不可怕，唯一可怕的是恰好完全平方！