素性测试与素数生成：Miller-Rabin、BPSW 与 OpenSSL/FIPS 186-5

关于素性测试，流传着三种说法：Miller-Rabin 要跑 64 轮才能把误判率压到 \(2^{-128}\)；竞赛里常用的 7 个底 \(\{2,3,5,7,11,13,37\}\) 覆盖全部 64 位整数；OpenSSL 生成 RSA 素数时用的是 Baillie-PSW 加增量搜索。三句话都不对。OpenSSL 3.6.2 生成 2048 位 RSA 密钥时，每个 1024 位素数只做 5 轮 Miller-Rabin，依据是 FIPS 186-5 的平均情况误差分析；那 7 个底在 \(2^{64}\) 以下有 52 个合数漏网，最小的是 4498414682539051；OpenSSL 的 BN_check_prime 里没有 Lucas 测试，素数生成循环在一次 Miller-Rabin 失败后会重新抽一个随机数，而不是继续往后找。

本文先用实验把 Miller-Rabin 的两种误差讲清楚：对任意合数成立的最坏情况 \(4^{-t}\)，和只对随机候选成立、小得多的平均情况界；再用 \(2^{64}\) 以下全部 2-伪素数表验证确定性底数集合；然后梳理 BPSW、AKS 的谱系；最后逐行对照 OpenSSL 3.6.2 源码和 FIPS 186-5，说明工业级素数生成到底做了什么。

复现程序都在文章目录的 reproduce/ 下，编译运行命令和外部数据的下载地址见其中的 README.md：

• 环境：Intel Core i9-12900K，WSL2 内核 6.6.87.2，GCC 16.1.1 20260430，GMP 6.3.0，Python 3.14.5；OpenSSL 源码与命令行均为 3.6.2（2026-04-07）。
• C 程序统一用 gcc -O2 -Wall -Wextra 编译，并在 -fsanitize=address,undefined 下运行过（规模见 README.md）。
• 所有指标都是计数（候选数、模幂次数、漏判个数），与时钟无关；随机程序用固定种子 20260923，重复运行输出相同。

一、试除与小素数筛：能排除多少候选

试除法（trial division）检查 \(2\) 到 \(\lfloor\sqrt n\rfloor\) 之间有没有因子。若 \(n=ab\) 且 \(1<a\le b\)，则 \(a\le\sqrt n\)，所以检查到 \(\sqrt n\) 就够了。代价是 \(O(\sqrt n)\) 次除法；对 \(k\) 位整数这是 \(2^{k/2}\) 量级，对输入长度是指数的。1024 位整数要试到 \(2^{512}\)，完全不可行。

试除在大数素性测试里的作用是预过滤：先用前几十到几千个小素数排除掉大部分合数，再对剩下的候选跑昂贵的模幂。能排除多少，有精确答案。随机奇数不被任何奇素数 \(p\le B\) 整除的概率是

\[ \prod_{3\le p\le B}\left(1-\frac1p\right)\approx\frac{2e^{-\gamma}}{\ln B}, \]

右边来自 Mertens 第三定理 \(\prod_{p\le B}(1-1/p)\sim e^{-\gamma}/\ln B\)，乘 2 是因为奇数已经排除了因子 2，\(\gamma\approx0.5772\) 是 Euler 常数。reproduce/prime_gen_sim.c 按乘积精确计算：

“\(k\) 个素数”包含 2，所以奇数实际只用 3 到 \(B\) 试除，与 OpenSSL 从 primes[1] 开始循环的写法一致。对 \(10^6\) 个随机 1024 位奇数做实测（用它们与 3 到 719 之积的最大公约数判断），存活 169,993 个，即 0.1700，与精确值 0.1697 一致。

表里能看到收益递减：从 128 个素数加到 2048 个，存活率只从 0.1697 降到 0.1146，筛的工作量却涨了 16 倍。存活率按 \(1/\ln B\) 衰减，再多的小素数也只能把候选压到一成左右，剩下的必须交给概率测试。

二、Fermat 测试与 Carmichael 数

Fermat 小定理：\(p\) 为素数且 \(\gcd(a,p)=1\) 时，\(a^{p-1}\equiv1\pmod p\)。逆否命题给出一个测试：若 \(a^{n-1}\not\equiv1\pmod n\)，则 \(n\) 必为合数，\(a\) 称为 Fermat 证人（witness）；若合数 \(n\) 满足 \(a^{n-1}\equiv1\)，\(a\) 称为 Fermat 说谎者（liar），\(n\) 称为以 \(a\) 为底的伪素数（pseudoprime）。

问题在于 Carmichael 数（Carmichael number）：对所有与 \(n\) 互素的 \(a\) 都有 \(a^{n-1}\equiv1\pmod n\) 的合数。Korselt（1899）给出了判据：合数 \(n\) 是 Carmichael 数，当且仅当 \(n\) 无平方因子，且对每个素因子 \(p\) 都有 \((p-1)\mid(n-1)\)。Alford、Granville、Pomerance（Annals of Mathematics, 1994）证明了 Carmichael 数有无穷多个。

reproduce/liars.c 对 20000 以下的每个奇合数枚举全部底数，找到 9 个 Carmichael 数：561、1105、1729、2465、2821、6601、8911、10585、15841，与 Korselt 判据逐个吻合。对它们，只有与 \(n\) 不互素的底才能揭穿合数身份。如果 Carmichael 数的素因子都很大，这样的底极少，多跑几轮 Fermat 测试也无济于事。

三、Miller-Rabin：强伪素数与 1/4 界

从 1 的平方根到强可能素数

Miller-Rabin 在 Fermat 测试上加了一条约束：模素数 \(p\) 时，\(x^2\equiv1\) 只有 \(x\equiv\pm1\) 两个解，因为 \(p\mid(x-1)(x+1)\) 意味着 \(p\) 整除其中一个因子。

把 \(n-1\) 写成 \(2^s d\)（\(d\) 为奇数）。

命题：若 \(n\) 为奇素数，则对任意 \(1\le a\le n-1\)，要么 \(a^d\equiv1\pmod n\)，要么存在 \(0\le r<s\) 使 \(a^{2^r d}\equiv-1\pmod n\)。

证明：考察序列 \(a^d, a^{2d}, \dots, a^{2^s d}=a^{n-1}\)，每一项是前一项的平方，最后一项由 Fermat 小定理等于 1。若 \(a^d\equiv1\)，命题成立；否则设 \(a^{2^r d}\) 是序列中第一个等于 1 的项，\(r\ge1\)，则 \(x=a^{2^{r-1}d}\) 满足 \(x^2\equiv1\) 且 \(x\not\equiv1\)，只能 \(x\equiv-1\)。\(\square\)

满足上述条件的奇合数称为以 \(a\) 为底的强伪素数（strong pseudoprime），\(a\) 称为强说谎者（strong liar）。测试过程就是算出 \(x_0=a^d\bmod n\)，再反复平方至多 \(s-1\) 次：遇到 \(n-1\) 判通过；遇到 1 而前一项不是 \(\pm1\)，或者平方到头也没出现 \(n-1\)，判合数。

图中第二行是 Fermat 测试漏掉的情形：561 对底 2 满足 \(2^{560}\equiv1\)，但平方链在 67 处出现了 1 的非平凡平方根。非平凡平方根还直接给出因子：\(67^2-1=(67-1)(67+1)\equiv0\)，于是 \(\gcd(66,561)=33\)、\(\gcd(68,561)=17\)。第三行说明强伪素数确实存在：\(2^{11}=2048\equiv1\pmod{2047}\)，而 \(1023=11\times93\)，所以 \(2^{1023}\equiv1\)，最小的以 2 为底的强伪素数就是 \(2047=23\times89\)。

Monier–Rabin 界：至多四分之一的底会说谎

定理（Monier 1980；Rabin 1980）：奇合数 \(n>9\) 的强说谎者个数 \(S(n)\) 满足 \(S(n)\le\varphi(n)/4\)，其中 \(\varphi\) 是 Euler 函数。

Monier 给出了 \(S(n)\) 的精确公式。设 \(n=\prod_{i=1}^{\omega}p_i^{e_i}\)，\(p_i-1=2^{s_i}d_i\)，\(\nu=\min_i s_i\)，\(d\) 为 \(n-1\) 的奇数部分，则

\[ S(n)=\left(1+\frac{2^{\omega\nu}-1}{2^{\omega}-1}\right)\prod_{i=1}^{\omega}\gcd(d,d_i). \]

以 \(561=3\cdot11\cdot17\) 为例：\(\omega=3\)，\(s_i=1,1,4\)，所以 \(\nu=1\)；\(d=35\)，\(d_i=1,5,1\)，\(\prod\gcd=5\)。代入得 \(S(561)=2\times5=10\)，而 \(\varphi(561)=320\)，即只有 3.1% 的底会说谎；作为对比，561 的 Fermat 说谎者是全部 320 个与它互素的底。reproduce/monier_check.py 用这个公式核对了 liars.c 暴力枚举的全部 7738 个奇合数，没有一个不符。定理的证明就是按 \(\omega=1\)（素数幂）、\(\omega=2\)、\(\omega\ge3\) 分情况估计这个公式。

实测结果（liars 20000）：

• \(n=9\) 时 \(S/\varphi=2/6=1/3\)，这是定理要求 \(n>9\) 的原因（说谎的底是 1 和 8）；
• \(n>9\) 时比值最大恰为 \(1/4\)，取等号的有 15、91、703、1891、8911、12403 六个数；
• 95.0% 的奇合数（7738 个中的 7353 个）强说谎者不到 \((n-1)/100\)。

所以 \(1/4\) 是一个很少被取到的最坏情况界。对固定的合数 \(n\)，独立随机选 \(t\) 个底全部说谎的概率不超过 \(4^{-t}\)，这对任何输入都成立，包括对手精心构造的合数。

谱系

• Miller（1976）：在扩展黎曼猜想（ERH）下给出确定性多项式时间判定。
• Solovay–Strassen（SIAM J. Comput., 1977）：基于 Euler 判别法和 Jacobi 符号的随机测试，每轮误判上界为 \(1/2\)。
• Rabin（1980）：去掉 ERH，改为随机选底，得到无条件的 \(1/4\) 界。
• Monier（1980）：独立给出同样的界，并得到上面的精确计数公式。

今天的”Miller-Rabin”指的就是 Rabin 的随机化版本。

四、平均情况：为什么 FIPS 只要 4 到 5 轮

\(4^{-t}\) 回答的是”一个给定的合数通过 \(t\) 轮的概率”。生成素数时关心的是另一个量：一个随机的 \(k\) 位奇数通过了 \(t\) 轮，它仍是合数的概率 \(p_{k,t}\)。FIPS 186-5 附录 C.1 专门用 WARNING 段区分这两者：前者对验证他人提供的数重要，后者对生成者重要。

两者差距很大，原因有二：随机合数的强说谎者比例通常远低于 1/4（上一节 95% 的合数不到 1%）；而且通过测试的数里素数占了绝大多数。Damgård、Landrock、Pomerance（Math. Comp., 1993，下称 DLP）给出了显式上界，FIPS 186-5 附录 C.1 把它写成公式 (2)：对 \(3\le M\le 2\sqrt{k-1}-1\)，

\[ p_{k,t}\le 2.00743\ln2\cdot k\cdot2^{-k}\left(2^{k-2-Mt}+\frac{8(\pi^2-6)}{3}\,2^{k-2}\sum_{m=3}^{M}2^{m-(m-1)t}\sum_{j=2}^{m}\frac{1}{2^{\,j+(k-1)/j}}\right), \]

对 \(M\) 取最小值，再找满足目标概率的最小 \(t\)。reproduce/fips_c1.py 用 log-sum-exp 形式实现这个算法，重算的轮数与 FIPS 186-5 Table B.1 全部一致：

1024 位时，\(t=1,2,3,4,5\) 对应的 \(\log_2 p_{1024,t}\) 分别是 \(-42.9\)、\(-72.8\)、\(-93.1\)、\(-109.8\)、\(-124.2\)，而最坏情况界只有 \(-2\)、\(-4\)、\(-6\)、\(-8\)、\(-10\)。位数越大，平均情况越好：随机大整数的强说谎者比例更低。FIPS 186-5 还指出，在同样的随机候选假设下，\(k\ge51\) 时可以证明 \(p_{k,t}\le4^{-t}\)。所以取 \(t=-\log_2(p_{\text{target}})/2\) 轮，对生成和验证都够用，只是保守。

用同一程序还能复核 OpenSSL 文档（doc/man3/BN_generate_prime.pod）的说法。OpenSSL 生成素数时至少做 64 轮，文档声称 512、1024、2048 位素数的误判率分别为 \(2^{-287}\)、\(2^{-435}\)、\(2^{-648}\)，更大的素数（128 轮）低于 \(2^{-882}\)。算出来是 \(2^{-287.5}\)、\(2^{-435.9}\)、\(2^{-648.2}\)，2049 位 128 轮为 \(2^{-883.2}\)，全部吻合。

这个平均情况在模拟里也看得到。第八节的生成实验一共测试了约 29 万个进入 Miller-Rabin 的 1024 位合数候选，没有一个通过哪怕第 1 轮。

五、64 位整数的确定性判定

最小强伪素数 \(\psi_k\) 与底数集合

对有界的 \(n\)，可以选一组固定底，使测试变成确定性的。记 \(\psi_k\) 为以前 \(k\) 个素数为底都通过的最小强伪素数，\(\psi_k\) 以下用前 \(k\) 个素数做底就不会出错。

检验”某组底在 \(2^{64}\) 以下确定”有一个干净的办法：以 2 为底的强伪素数都是以 2 为底的 Fermat 伪素数，而 Jan Feitsma 计算、William Galway 整理发布了 \(2^{64}\) 以下全部 2-Fermat 伪素数的列表。只要底数集合含 2，逐个检查列表里的数有没有通过全部底即可。reproduce/psp64_verify.c 读入整张表：共 118,968,378 个 2-伪素数，其中 31,894,014 个是 2-强伪素数。算出的 \(\psi_k\) 与文献一致：

最后两行超出 \(2^{64}\approx1.84\times10^{19}\)，程序只能验证”\(2^{64}\) 以下不存在”；它们的精确值引自 Sorenson 与 Webster。常见的错误是把 \(3.3\times10^{24}\) 写成 12 个底的适用范围，那其实是 13 个底的界；12 个底的界是 \(\psi_{12}\approx3.19\times10^{23}\)，依然远大于 \(2^{64}\)。

几组常见底数集合在 \(2^{64}\) 以下的漏判数：

\(\{2,7,61\}\) 在 4759123141 以下是确定的，这个数大于 \(2^{32}\)，所以它适合 32 位整数。\(\{2,3,5,7,11,13,37\}\) 的最小反例 \(4498414682539051=46411\times232051\times417691\) 通过了全部 7 个底，它不能用于 64 位。Jim Sinclair 在 2011 年找到的 7 底集合（发布在 miller-rabin.appspot.com）经整表验证没有漏判，比前 12 个素数少 5 次模幂。

实现要点

reproduce/mr64.h 的核心是 20 来行：

typedef unsigned __int128 u128;

static inline uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m)
{
    return (uint64_t)((u128)a * b % m);
}

static inline bool sprp(uint64_t n, uint64_t a)   /* n odd, n > 2 */
{
    uint64_t d = n - 1;
    int s = __builtin_ctzll(d);
    d >>= s;

    a %= n;
    if (a == 0)
        return true;

    uint64_t x = powmod(a, d, n);
    if (x == 1 || x == n - 1)
        return true;
    for (int r = 1; r < s; r++) {
        x = mulmod(x, x, n);
        if (x == n - 1)
            return true;
        if (x == 1)
            return false;
    }
    return false;
}

有两处容易写错，reproduce/mr64_test.c 各自构造了反例：

• 乘法溢出。两个接近 \(2^{64}\) 的数相乘需要 128 位中间结果。把 mulmod 换成 64 位的 a * b % m，\((2^{64}-20059,\,2^{64}-59]\) 里的 447 个素数全部被判为合数。GCC 和 Clang 在 64 位目标上提供 unsigned __int128；没有它时要用 Montgomery 乘法（见下一篇）。
• 底数是 \(n\) 的倍数。Sinclair 的底很大，素数 \(n\) 可能整除某个底：\(28178=2\cdot73\cdot193\)，\(450775=5^2\cdot13\cdot19\cdot73\)，\(9780504=2^3\cdot3\cdot407521\)，\(1795265022=2\cdot3\cdot299210837\)。此时 \(a\bmod n=0\)，\(0^d=0\) 既不是 1 也不是 \(n-1\)。不加 a == 0 判断时，73、193、407521、299210837 这四个素数都会被判为合数。另一种办法是先用这些小素数试除。

交叉测试结果：

• \(10^8\) 以下逐个与埃氏筛比对，前 12 素数版、Sinclair 版、BPSW 版全部一致，\(\pi(10^8)=5761455\)；
• 200 万个随机 64 位奇数（其中 92,092 个是素数），加上 \(2^{32}\)、\(2^{63}\)、\(2^{64}\) 附近三个各 20 万个数的窗口，与 GMP 6.3.0 的 mpz_probab_prime_p 全部一致；
• 表中的各个 \(\psi_k\) 恰好通过前 \(k\) 个底；
• 最大的 64 位素数是 \(2^{64}-59=18446744073709551557\)。

在 GRH 下的确定性

Miller（JCSS, 1976）证明：若扩展黎曼猜想成立，最小的 Miller-Rabin 证人是 \(O((\ln n)^2)\)，由此得到 \(O((\log n)^4)\) 的确定性判定，但他没有给出常数。Bach（Math. Comp., 1990）给出显式界：在 ERH 下，每个奇合数都有小于 \(2(\ln n)^2\) 的证人。对 \(n<2^{64}\)，这个界约为 3936 个底，远不如上面 7 个底的整表验证结果。它的意义是理论上的：如果 ERH 成立，素性判定就有简单的确定性多项式算法。

六、Baillie-PSW：两类伪素数几乎不相交

组合方式

Pomerance、Selfridge、Wagstaff（Math. Comp. 35(151), 1980）与 Baillie、Wagstaff（Math. Comp. 35(152), 1980）提出把两种结构不同的测试组合起来，后人称为 Baillie-PSW（BPSW）：

1. 以 2 为底的强可能素数测试；
2. 强 Lucas 可能素数测试（strong Lucas probable prime test）。

Lucas 测试用的是 Lucas 序列 \(U_k,V_k\)，它们由 \(x^2-Px+Q\) 的两个根 \(\alpha,\beta\) 定义：\(U_k=(\alpha^k-\beta^k)/(\alpha-\beta)\)，\(V_k=\alpha^k+\beta^k\)。取判别式 \(D=P^2-4Q\) 满足 Jacobi 符号 \(\left(\frac{D}{n}\right)=-1\)，此时对素数 \(n\) 有 \(n\mid U_{n+1}\)。把 \(n+1\) 写成 \(2^s d\)，强 Lucas 测试要求 \(U_d\equiv0\)，或存在 \(0\le r<s\) 使 \(V_{2^r d}\equiv0\pmod n\)。参数按 Selfridge 的方法 A 选取：\(D\) 依次取 \(5,-7,9,-11,\dots\) 中第一个满足 \(\left(\frac{D}{n}\right)=-1\) 的值，\(P=1\)，\(Q=(1-D)/4\)。完全平方数找不到这样的 \(D\)，要先排除。

直观上，Miller-Rabin 检验的是乘法群 \((\mathbb Z/n\mathbb Z)^*\) 里与 \(n-1\) 有关的结构，Lucas 测试检验的是二次扩张里与 \(n+1\) 有关的结构。Baillie 与 Wagstaff（1980，第 1404–1405 页）观察到，2-伪素数和 Lucas 伪素数对小模数分别倾向于落在 \(+1\) 和 \(-1\) 的剩余类里，两类伪素数因此很难重合。

mr64.h 实现了这个测试。强 Lucas 伪素数单独来看并不少：100000 以下有 12 个（5459、5777、10877、16109、18971、22499、24569、25199、40309、58519、75077、97439，与 SymPy 的 is_strong_lucas_prp 一致），\(10^8\) 以下有 505 个。但 \(2^{64}\) 以下 31,894,014 个 2-强伪素数中，没有一个同时通过强 Lucas 测试，即 BPSW 在 \(2^{64}\) 以下没有反例。Baillie、Fiori、Wagstaff（2021）报告了更强的结论：118,968,378 个 2-伪素数中没有一个是 Lucas 伪素数。他们用的参数方法 A* 在 \(Q=-1\) 时把 \(P,Q\) 都改成 5，并证明方法 A 与 A* 产生的 Lucas 伪素数和强 Lucas 伪素数完全相同，所以两者的结果可以直接对照。

反例：没找到，但预计存在

BPSW 的作者为反例设立的 620 美元悬赏，40 年后仍无人领取（Baillie–Fiori–Wagstaff 2021）。另一方面，Pomerance 在 1984 年的短文 “Are there counterexamples to the Baillie-PSW primality test?” 中改造了 Erdős 构造 Carmichael 数的启发式论证：反例应当有无穷多个，只是没有构造方法，也不知道最小的有多大。这是素性测试里最有名的开放问题之一，第十节再谈。

谁在用

• GMP 6.3.0 手册写明，mpz_probab_prime_p(n, reps) 先做若干试除，再做一次 BPSW，然后做 \(\text{reps}-24\) 轮 Miller-Rabin。
• FIPS 186-5 附录 B.3 允许在 Miller-Rabin 之后加一次 Lucas 测试，理由是”没有已知的合数能通过这一组合”，但不强制。
• OpenSSL 3.6.2 的 BN_check_prime 不做 Lucas 测试，只做试除和随机底 Miller-Rabin（见第八节）。

七、确定性证明：AKS、ECPP 与 Lucas–Lehmer

AKS

Agrawal、Kayal、Saxena 在 2002 年发布预印本 “PRIMES is in P”，正式论文发表于 Annals of Mathematics 160(2), 2004，给出第一个无条件的确定性多项式时间素性判定，并获得 2006 年 Gödel 奖和 Fulkerson 奖。出发点是多项式恒等式：对与 \(n\) 互素的 \(a\)，\(n\) 为素数当且仅当

\[ (X+a)^n\equiv X^n+a\pmod n. \]

直接验证要展开 \(n+1\) 个系数。AKS 改为在模 \(X^r-1\) 和 \(n\) 的意义下，对 \(O(\sqrt r\log n)\) 个 \(a\) 验证，其中 \(r=O(\log^5 n)\) 量级，从而得到多项式时间。

复杂度的说法常被写错：

• 预印本证明的是 \(\tilde O(\log^{12}n)\)；
• 正式论文证明 \(\tilde O(\log^{21/2}n)\)；借助一个常数不可计算的解析数论结果，可以改进到 \(\tilde O(\log^{15/2}n)\)；
• 在一个关于 Sophie Germain 素数密度的猜想下是 \(\tilde O(\log^6 n)\)；
• Lenstra 与 Pomerance（J. Eur. Math. Soc. 21(4), 2019）用 Gaussian 周期构造了一个变体，无条件证明了 \((\log n)^6(2+\log\log n)^c\)。

AKS 的价值在于理论，密码学库生成或验证素数时并不用它。

ECPP

需要可验证的素性证书时用的是椭圆曲线素性证明（elliptic curve primality proving，ECPP，Atkin–Morain, Math. Comp., 1993）。它输出一串证书，别人可以比生成快得多地验证；运行时间的分析依赖启发式假设。本文不展开。

Lucas–Lehmer

Mersenne 数 \(M_p=2^p-1\)（\(p\) 为奇素数）有专用判定：令 \(S_0=4\)，\(S_i=S_{i-1}^2-2\)，则 \(M_p\) 为素数当且仅当 \(S_{p-2}\equiv0\pmod{M_p}\)。

它快的原因常被说错。Lucas–Lehmer 需要 \(p-2\) 次模 \(M_p\) 平方；对 \(M_p\) 做一次 Fermat 或 Miller-Rabin 测试，模幂的指数有 \(p\) 位，同样约需 \(p\) 次平方。两者的平方次数相当，Lucas–Lehmer 的优势是确定性，而不是更少的运算。真正让 Mersenne 数好算的是：

• 模 \(2^p-1\) 的约简只需移位加法：\(x\equiv(x\bmod2^p)+\lfloor x/2^p\rfloor\)；
• 大数平方可以用 FFT 乘法（见上一篇）。

GIMPS 当前的流程也说明了这一点。客户端先做 Fermat 可能素数（PRP）测试，出现候选后再做确认。已知最大的 Mersenne 素数 \(M_{136279841}\) 有 41,024,320 位十进制数字：2024 年 10 月 11 日，一次 PRP 测试（NVIDIA A100 上）报告它可能是素数；10 月 12 日，Lucas–Lehmer 测试（H100 上）确认，官方发现日期记为 10 月 12 日（mersenne.org 公告）。

八、OpenSSL 3.6.2 怎么判素、怎么生成素数

以下源码引自 OpenSSL 3.6.2（tag openssl-3.6.2），路径相对于源码根目录。

BN_check_prime：按最坏情况设计

BN_is_prime_ex 与 BN_is_prime_fasttest_ex 在 3.0 被标为弃用（CHANGES.md，Kurt Roeckx），取代它们的是 3.0 新增的 BN_check_prime(p, ctx, cb)，它没有轮数参数。crypto/bn/bn_prime.c 中决定规模的两个函数是：

static int calc_trial_divisions(int bits)
{
    if (bits <= 512)
        return 64;
    else if (bits <= 1024)
        return 128;
    else if (bits <= 2048)
        return 384;
    else if (bits <= 4096)
        return 1024;
    return NUMPRIMES;           /* 2048 */
}

/*
 * Use a minimum of 64 rounds of Miller-Rabin, which should give a false
 * positive rate of 2^-128. If the size of the prime is larger than 2048
 * the user probably wants a higher security level than 128, so switch
 * to 128 rounds giving a false positive rate of 2^-256.
 */
static int bn_mr_min_checks(int bits)
{
    if (bits > 2048)
        return 128;
    return 64;
}

BN_check_prime 调用 ossl_bn_check_prime(p, 0, ctx, 1, cb)，后者把轮数钳到至少 bn_mr_min_checks()，然后进入 bn_is_prime_int()：先用 primes[1] 到 primes[K-1] 试除（\(K\) 取自 calc_trial_divisions），再调用 ossl_bn_miller_rabin_is_prime()，每轮用 BN_priv_rand_range_ex 在 \([2, w-2]\) 中取随机底。弃用的两个函数也走同一个 ossl_bn_check_prime，传入更小的轮数同样会被钳到 64 或 128。注释里的 \(2^{-128}\) 就是 \(4^{-64}\)，是最坏情况界。检查外来的数（例如对方提供的 DH 参数）需要的正是这一种。

include/openssl/bn.h 里还保留着弃用的宏 BN_prime_checks_for_size(b)：3747 位以上 3 轮，1345 位以上 4 轮，476 位以上 5 轮，依次到 55 位以上 27 轮、更小的 34 轮。注释说这张表按 FIPS 186-4 附录 F.1 的算法生成，是平均情况的轮数，3.x 的判素函数已经不再使用它。

BN_generate_prime_ex2：筛，失败就重抽

生成循环（crypto/bn/bn_prime.c，删去了 safe prime 与 add/rem 分支）：

    int checks = bn_mr_min_checks(bits);
    /* ... */
loop:
    /* make a random number and set the top and bottom bits */
    if (!probable_prime(ret, bits, safe, mods, ctx))
        goto err;
    if (!BN_GENCB_call(cb, 0, c1++))
        goto err;
    i = bn_is_prime_int(ret, checks, ctx, 0, cb);   /* no trial division */
    if (i == -1)
        goto err;
    if (i == 0)
        goto loop;

probable_prime() 用 BN_priv_rand_ex(rnd, bits, BN_RAND_TOP_TWO, BN_RAND_BOTTOM_ODD, ...) 取一个最高两位为 1 的随机奇数，算出它对每个小素数的余数 mods[i]，然后让 delta 从 0 起每次加 2，直到所有 (mods[i] + delta) % primes[i] 都不为 0，返回 rnd + delta。初始化 mods[i] 时，每个小素数要做一次大数除以机器字；之后每一步只做机器字运算。Miller-Rabin 失败后，goto loop 会重新调用 probable_prime()，也就是重新抽一个随机数。增量搜索只发生在筛的内部，不跨越 Miller-Rabin 失败。文件里的注释称这种快速筛法来自 PGP 中 Philip Zimmermann 的实现。

图里的例子是真实数据：1000031 的最小素因子 41 超出了筛的范围，只能靠 Miller-Rabin 拒绝，它的 1000030 个底中只有 50 个是强说谎者。

reproduce/prime_gen_sim.c 按这个流程（GMP 大数，随机底，固定种子）生成 1024 位素数，统计每得到一个素数的开销：

除第一行外每行生成 1000 个素数，每个素数做 64 轮。估计值的来历：\(k\) 位随机奇数是素数的概率约为 \(2/\ln N\)，筛后存活的比例是 \(s\)（第一节的存活率），取 \(N\approx0.875\cdot2^{1024}\)（最高两位为 1 时的平均值）。实测与估计相差最多 6%，约为两倍抽样标准误（每组候选数近似几何分布，1000 个样本的标准误约为均值的 3%）。常见的”1024 位素数平均要试约 710 个候选”把偶数也算进去了；只取奇数时约为 355 个。

模幂次数揭示了轮数的真实代价。几乎所有合数候选在第 1 轮就被拒绝，所以每个素数的开销约为（合数候选数）+（轮数）。同样用 128 个素数筛：5 轮时每个素数 64.0 次模幂，64 轮时 125.6 次。64 轮比 5 轮贵约一倍，而不是 13 倍。

RSA 密钥生成走另一条路

openssl genrsa 2048 并不调用 BN_generate_prime_ex2。crypto/rsa/rsa_gen.c 的 rsa_keygen() 在非 FIPS 构建中这样分派：

    if (primes == 2
        && bits >= 2048
        && (e_value == NULL || BN_num_bits(e_value) > 16))
        ok = ossl_rsa_sp800_56b_generate_key(rsa, bits, e_value, cb);
    else
        ok = rsa_multiprime_keygen(rsa, bits, primes, e_value, cb);

\(65537=2^{16}+1\) 有 17 位，所以常见的 2048 位、\(e=65537\) 密钥走 SP 800-56B 路径（FIPS 构建中则总是走这条路径）。CHANGES.md 在 3.0 条目中记载，双素数 RSA 的默认生成方法改为 FIPS 186-4 B.3.6（Shane Lontis）。

flowchart TD
    A["rsa_keygen()"] --> B{"2 primes, bits >= 2048, e longer than 16 bits?"}
    B -- "yes, e.g. 2048-bit key with e = 65537" --> C["ossl_rsa_sp800_56b_generate_key()"]
    B -- "no" --> D["rsa_multiprime_keygen()"]
    D --> E["BN_generate_prime_ex2() per prime: sieve, 64 MR rounds (128 above 2048 bits)"]
    C --> F["aux primes p1, p2: random odd start of 141 / 171 / 201 bits, step +2"]
    F --> G["each step: trial division, then 38 / 41 / 44 MR rounds"]
    G --> H["Y = X + ((R - X) mod 2*r1*r2), then step Y by 2*r1*r2, at most 20 * nlen/2 steps"]
    H --> I["each step: gcd(Y - 1, e) = 1, trial division, 5 MR rounds (4 if nlen >= 3072)"]
    I --> J["repeat for q; require |Xp - Xq| and |p - q| > 2^(nlen/2 - 100)"]

crypto/bn/bn_rsa_fips186_4.c 中的轮数直接写着出处：

/*
 * Refer to FIPS 186-5 Table B.1 for minimum rounds of Miller Rabin
 * required for generation of RSA primes (p and q)
 */
static int bn_rsa_fips186_5_prime_MR_rounds(int nbits)
{
    if (nbits >= 3072)
        return 4;
    if (nbits >= 2048)
        return 5;
    return 0; /* Error */
}

辅助素数的轮数由 bn_rsa_fips186_5_aux_prime_MR_rounds() 给出，nlen 为 2048、3072、4096 位时分别是 38、41、44 轮，最小位长分别是 141、171、201 位。候选检测调用的是 ossl_bn_check_generated_prime()：它总做试除，轮数照原样使用，不钳到 64。注释还说明，候选上限 imax 从 FIPS 186-4 的 \(5\cdot\text{nlen}/2\) 提高到 FIPS 186-5 附录 B.9 的 \(20\cdot\text{nlen}/2\)。据 Allen Roginsky 的分析，旧上限约有两百万分之一的失败概率。

genrsa -verbose 的进度字符可以直接印证这些轮数。apps/lib/s_cb.c 的 progress_cb 按回调事件 0、1、2、3 分别输出 .、+、*、换行。在本机运行三次 openssl genrsa -verbose -out /dev/null 2048，每次输出中长串 + 的长度都是 39、39、6、39、39、6。下面是其中一次的开头（截断）：

Generating RSA key with 2048 bits
............+...+........+++++++++++++++++++++++++++++++++++++++*.+.+..+.......+.....+.+...

对照源码：

• 每个 . 是一个新候选（事件 0）；
• 单个 + 是候选通过了试除（bn_is_prime_int 的事件 1，参数 \(-1\)），但在第 1 轮 Miller-Rabin 失败；
• 辅助素数的 39 个 + 是 1 次试除加 38 轮，之后的 * 是事件 2；
• \(p\)、\(q\) 的 6 个 + 是 1 次试除加 5 轮，之后的换行是事件 3。

所以 OpenSSL 3.6.2 生成 2048 位 RSA 密钥时，\(p\) 和 \(q\) 各只做 5 轮 Miller-Rabin，没有 Lucas 测试。

九、FIPS 186-5 对 RSA 素数的要求

FIPS 186-5（2023）附录 A.1.1 对 IFC（整数分解类）密钥规定：

• \(e\) 为奇数，\(2^{16}<e<2^{256}\)，在生成 \(p,q\) 之前选定；
• \(p,q\) 随机选取，\(\sqrt2\cdot2^{\text{nlen}/2-1}\le p,q\le2^{\text{nlen}/2}-1\)，且 \(p-1\)、\(q-1\) 与 \(e\) 互素；
• \(|p-q|>2^{\text{nlen}/2-100}\)；
• \(2^{\text{nlen}/2}<d<\operatorname{LCM}(p-1,q-1)\)，极少数不满足时重新生成。

\(p,q\) 的生成方法有两大类：随机素数（A.1.2 可证素数、A.1.3 可能素数），以及带条件的素数（A.1.4 至 A.1.6）。带条件的素数要求 \(p-1\)、\(p+1\)、\(q-1\)、\(q+1\) 分别有辅助素因子 \(p_1,p_2,q_1,q_2\)，其最小长度见 Table A.1：nlen 在 2048 到 3071 之间时大于 140 位，3072 到 4095 之间时大于 170 位，4096 及以上时大于 200 位。OpenSSL 走的是 A.1.6（全部为可能素数）。

A.1.3 的过程值得对照 OpenSSL：

• nlen 小于 2048 直接返回 FAILURE；
• 每个候选都是从 DRBG 新取的 nlen/2 位串，可以把最高两位置 1，低位调成奇数，小于 \(\sqrt2\cdot2^{\text{nlen}/2-1}\) 就重取；
• 满足 \(\gcd(p-1,e)=1\) 的候选按 Table B.1 的轮数测试；
• 尝试 \(5\cdot\text{nlen}\) 次仍未找到 \(p\)（\(q\) 为 \(10\cdot\text{nlen}\) 次）即失败。

也就是说，标准的”随机可能素数”方法每个候选都重新抽样，不做增量搜索。

附录 B.3 规定：

• 试除可选，上限 \(L\) 取 \(10^3\) 到 \(10^5\) 之间；
• Miller-Rabin 的最少轮数见 Table B.1（第四节已全部重算）；
• 可以追加一次 Lucas 测试，辅助素数不需要。

与 FIPS 186-4（2013）相比，186-4 Table C.2、C.3 还列有 512 位 \(p,q\) 的行（例如 \(2^{-80}\) 目标下 5 轮，辅助素数 28 轮），186-5 的生成过程则直接拒绝 2048 位以下的 nlen。

附录 C.2 解释了辅助素数的用意：选择它们的长度，是为了让”资源可观但不无限”的攻击者用 Lenstra 的椭圆曲线分解法（ECM）比用 Pollard \(p-1\)、Williams \(p+1\) 或各种循环攻击更划算；ECM 本身除了把 \(p,q\) 选大之外无法防御。

十、争论与开放问题

BPSW 会不会有反例

一方是 \(2^{64}\) 以下没有反例，四十多年来 620 美元的悬赏无人领取，第六节的实验对全部 31,894,014 个 2-强伪素数复核了这一点；另一方是 Pomerance（1984）的启发式论证：反例有无穷多个。两者并不矛盾：启发式预测的反例可能极大，远超任何穷举范围。Baillie、Fiori、Wagstaff（Math. Comp., 2021；预印本 arXiv:2006.14425）因此提出一个加强版，额外检查 \(V_{n+1}\equiv2Q\) 与 \(Q^{(n+1)/2}\equiv Q\cdot\left(\frac{Q}{n}\right)\pmod n\) 两个同余。他们在 \(10^{15}\) 以下只找到 5 个满足前一个同余的合数（Lucas-V 伪素数），而满足 \(U_{n+1}\equiv0\) 的 Lucas 伪素数有两百多万个。Baillie 与 Wagstaff 各自为加强版的第一个反例，或第一篇证明不存在反例的同行评审论文，悬赏 1000 美元；他们也在文中复述了 Pomerance 的论证，认为加强版同样有无穷多个反例。这个问题可以被证伪：只要有人给出一个具体反例。

对抗性输入：平均情况界不能用来验证别人的数

第四节的 \(p_{k,t}\) 只对随机候选成立。Albrecht、Massimo、Paterson、Somorovsky 在 “Prime and Prejudice”（ACM CCS 2018）中研究了输入由对手构造的情形：

• 他们构造的 2048 位合数能以 \(1/16\) 的概率通过 OpenSSL 当时默认配置下的素性测试，而文档宣称的误判率是 \(2^{-80}\)；
• 1024 位合数在 GMP 按推荐最少轮数配置时总能通过；
• Cryptlib、LibTomCrypt、JavaScript Big Number、WolfSSL 可以被总能通过的合数欺骗；
• 他们展示了这些漏洞如何被用来构造表面安全、实际离散对数容易的恶意 DH 参数，并建议使用 BPSW。

现在的 OpenSSL 3.6.2 在 BN_check_prime 里用至少 64 轮随机底（最坏情况 \(2^{-128}\)），GMP 6.3.0 在 Miller-Rabin 之前先做 BPSW。争论点在于：生成与验证是否应该用同一个函数、同一个轮数。OpenSSL 的选择是分开：验证按最坏情况，RSA 生成按 FIPS 的平均情况。

要不要”强素数”

Rivest 与 Silverman 在 “Are ‘Strong’ Primes Needed for RSA?”（IACR ePrint 2001/007）中认为，要求 \(p-1\)、\(p+1\) 有大素因子没有必要：对足够大的随机素数，Pollard \(p-1\) 和 Williams \(p+1\) 成功的概率本来就可以忽略，而 ECM 无法用这种条件防御。FIPS 186-5 保留了带辅助素数的生成方法，理由见附录 C.2；但它同时允许 A.1.2、A.1.3 的纯随机素数，也就是说标准本身没有把辅助素数当作必需。OpenSSL 选择了带条件的方法。

增量搜索的输出分布

从随机起点往后找第一个素数时，素数被选中的概率与它前面素数间隙的长度成正比，输出不是均匀分布。Brandt 与 Damgård（CRYPTO ’92）分析了这种增量搜索：给出了输出合数概率的显式上界；在 Hardy–Littlewood 素数 \(r\) 元组猜想下给出输出分布熵的下界，\(k\to\infty\) 时熵几乎达到最大值。允许重选起点的变体，他们只做了部分分析。OpenSSL 的两条路径都是混合流程：BN_generate_prime_ex2 只在筛内增量、Miller-Rabin 失败即重抽；FIPS 路径的辅助素数是从随机起点 +2 找到的第一个可能素数，\(p\) 则沿公差为 \(2r_1r_2\) 的等差数列搜索。它们的输出熵依赖一个未证明的猜想，本文没有测量。

其他开放问题

• 确定性底数集合还能不能再小。7 个底覆盖 \(2^{64}\) 是 2011 年的搜索结果；本文没有查到”6 个底不够”的证明，也没有查到反例。
• 无条件的确定性算法能否达到 Miller–Bach 在 ERH 下的 \(\tilde O(\log^4 n)\)。AKS 路线目前最好的无条件结果是 Lenstra–Pomerance 的 \(\tilde O(\log^6 n)\)。

十一、工程陷阱

选型可以归结为三种场景：

• 64 位以内：用确定性底数集合，前 12 个素数或 Sinclair 的 7 个底，乘法用 128 位。
• 判定外来的大数：用 BPSW，或者 BPSW 加若干轮随机底；不要依赖平均情况界。
• 生成密码学素数：候选必须来自可靠的 DRBG，筛掉小因子后按 FIPS 186-5 Table B.1 的轮数测试就足够。多加轮数只让每个素数的模幂次数从约 64 增加到约 126（1024 位，第八节），代价不大，但它防的是一个平均情况下已经小于 \(2^{-112}\) 的风险。

十二、参考资料

规范与文档

• NIST, FIPS 186-5, Digital Signature Standard (DSS), 2023：附录 A.1.1、A.1.3、A.1.6、Table A.1，附录 B.3、Table B.1、B.9，附录 C.1 公式 (1)(2) 与 WARNING 段，附录 C.2。
• NIST, FIPS 186-4, Digital Signature Standard (DSS), 2013：附录 C.3，Table C.2、C.3，附录 F.1。
• OpenSSL 3.6.2 手册 doc/man3/BN_generate_prime.pod（BN_check_prime、BN_generate_prime_ex2，默认轮数与误判率）。
• GMP 6.3.0 Manual, “Number Theoretic Functions”：mpz_probab_prime_p。

源码

• OpenSSL 3.6.2 crypto/bn/bn_prime.c：calc_trial_divisions()、bn_mr_min_checks()、BN_generate_prime_ex2()、ossl_bn_check_prime()、ossl_bn_check_generated_prime()、BN_check_prime()、bn_is_prime_int()、ossl_bn_miller_rabin_is_prime()、probable_prime()。
• OpenSSL 3.6.2 crypto/bn/bn_rsa_fips186_4.c：bn_rsa_fips186_5_aux_prime_MR_rounds()、bn_rsa_fips186_5_prime_MR_rounds()、bn_rsa_fips186_5_aux_prime_min_size()、ossl_bn_rsa_fips186_4_derive_prime()。
• OpenSSL 3.6.2 crypto/rsa/rsa_gen.c：rsa_keygen()；crypto/rsa/rsa_sp800_56b_gen.c：ossl_rsa_sp800_56b_generate_key() 等；include/openssl/bn.h：BN_prime_checks_for_size；apps/lib/s_cb.c：progress_cb()；CHANGES.md（3.0 条目）。

核心论文

• G. L. Miller, “Riemann’s Hypothesis and Tests for Primality”, Journal of Computer and System Sciences 13(3), 300–317, 1976.
• M. O. Rabin, “Probabilistic Algorithm for Testing Primality”, Journal of Number Theory 12(1), 128–138, 1980.
• L. Monier, “Evaluation and Comparison of Two Efficient Probabilistic Primality Testing Algorithms”, Theoretical Computer Science 12(1), 97–108, 1980.
• C. Pomerance, J. L. Selfridge, S. S. Wagstaff Jr., “The Pseudoprimes to \(25\cdot10^9\)”, Mathematics of Computation 35(151), 1003–1026, 1980.
• R. Baillie, S. S. Wagstaff Jr., “Lucas Pseudoprimes”, Mathematics of Computation 35(152), 1391–1417, 1980.
• I. Damgård, P. Landrock, C. Pomerance, “Average Case Error Estimates for the Strong Probable Prime Test”, Mathematics of Computation 61(203), 177–194, 1993.
• M. Agrawal, N. Kayal, N. Saxena, “PRIMES is in P”, Annals of Mathematics 160(2), 781–793, 2004.

其他论文

• A. Korselt, “Problème chinois”, L’Intermédiaire des Mathématiciens 6, 142–143, 1899.
• R. Solovay, V. Strassen, “A Fast Monte-Carlo Test for Primality”, SIAM Journal on Computing 6(1), 84–85, 1977.
• E. Bach, “Explicit Bounds for Primality Testing and Related Problems”, Mathematics of Computation 55(191), 355–380, 1990.
• G. Jaeschke, “On Strong Pseudoprimes to Several Bases”, Mathematics of Computation 61(204), 915–926, 1993.
• A. O. L. Atkin, F. Morain, “Elliptic Curves and Primality Proving”, Mathematics of Computation 61(203), 29–68, 1993.
• W. R. Alford, A. Granville, C. Pomerance, “There Are Infinitely Many Carmichael Numbers”, Annals of Mathematics 139(3), 703–722, 1994.
• Y. Jiang, Y. Deng, “Strong Pseudoprimes to the First Eight Prime Bases”, Mathematics of Computation 83, 2915–2924, 2014.
• J. Sorenson, J. Webster, “Strong Pseudoprimes to Twelve Prime Bases”, Mathematics of Computation 86(304), 985–1003, 2017.
• H. W. Lenstra Jr., C. Pomerance, “Primality Testing with Gaussian Periods”, Journal of the European Mathematical Society 21(4), 1229–1269, 2019.
• C. Pomerance, “Are There Counterexamples to the Baillie-PSW Primality Test?”, in H. W. Lenstra Jr., J. K. Lenstra, P. van Emde Boas (eds.), Dopo Le Parole, Amsterdam, 1984.
• R. Baillie, A. Fiori, S. S. Wagstaff Jr., “Strengthening the Baillie-PSW Primality Test”, Mathematics of Computation, 2021；预印本 arXiv:2006.14425。
• J. Brandt, I. Damgård, “On Generation of Probable Primes by Incremental Search”, CRYPTO ’92, LNCS 740, 358–370, 1993.
• M. R. Albrecht, J. Massimo, K. G. Paterson, J. Somorovsky, “Prime and Prejudice: Primality Testing Under Adversarial Conditions”, ACM CCS 2018；IACR ePrint 2018/749。
• R. L. Rivest, R. D. Silverman, “Are ‘Strong’ Primes Needed for RSA?”, IACR ePrint 2001/007.
• M. Nemec, M. Sýs, P. Švenda, D. Klinec, V. Matyáš, “The Return of Coppersmith’s Attack: Practical Factorization of Widely Used RSA Moduli”, ACM CCS 2017, 1631–1648.
• N. Heninger, Z. Durumeric, E. Wustrow, J. A. Halderman, “Mining Your Ps and Qs: Detection of Widespread Weak Keys in Network Devices”, USENIX Security 2012.

工程资料

• J. Feitsma, W. Galway, Tables of pseudoprimes and related data（\(2^{64}\) 以下全部 2-Fermat 伪素数），https://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html。
• J. Sinclair 的 7 底集合，见 https://miller-rabin.appspot.com/（2011-04-20 条目）。
• GIMPS, “Mersenne Prime Discovery - 2^136279841-1 is Prime!”, mersenne.org, 2024-10-21。

实验

• reproduce/mr64.h、reproduce/mr64_test.c：第五节的实现与交叉测试，第六节的强 Lucas 伪素数。
• reproduce/psp64_verify.c：第五节的 \(\psi_k\) 与底数集合表，第六节的 BPSW 复核。
• reproduce/liars.c、reproduce/monier_check.py、reproduce/plot_liars.py：第二、三节的说谎者统计与图。
• reproduce/fips_c1.py、reproduce/plot_error_bounds.py：第四节的 Table B.1 重算与图。
• reproduce/prime_gen_sim.c：第一节的试除表与第八节的生成统计。

系列导航： - 上一篇：快速傅里叶变换 - 下一篇：扩展欧几里得与模逆元

相关阅读： - RSA 从原理到攻击 - 有限域算术 - 椭圆曲线算术 - 随机化算法