素性测试与工业级素数生成

RSA 密钥生成的第一步是找到两个大素数 p 和 q。这听起来简单，但当数字达到 1024 位甚至 2048 位时，问题变得相当微妙：我们无法枚举所有因子，也无法用确定性算法在可接受的时间内完成判定。工业界的做法是用概率算法——特别是 Miller-Rabin 测试——在极高的置信度下完成素性判定。

本文从最朴素的试除法出发，逐步深入到 Fermat 小定理、Miller-Rabin 测试、Baillie-PSW 测试、AKS 确定性算法和 Lucas-Lehmer 测试，最终回到工程实践：OpenSSL 是如何生成素数的，RSA 密钥生成中有哪些需要注意的陷阱。

一、试除法及其局限

试除法是最直观的素性测试方法。给定一个正整数 n，只需检查 2 到 sqrt(n) 之间是否存在 n 的因子。

1.1 基本原理

如果 n 是合数，那么 n = a * b，其中 a 和 b 都大于 1。a 和 b 中至少有一个不超过 sqrt(n)，否则 a * b > n，产生矛盾。因此我们只需要检查到 sqrt(n) 即可。

1.2 C 实现

#include <stdbool.h>
#include <math.h>

bool is_prime_trial(unsigned long long n)
{
    if (n < 2)
        return false;
    if (n < 4)
        return true;
    if (n % 2 == 0 || n % 3 == 0)
        return false;

    /* 只需检查 6k +/- 1 形式的数 */
    for (unsigned long long i = 5; i * i <= n; i += 6) {
        if (n % i == 0 || n % (i + 2) == 0)
            return false;
    }
    return true;
}

1.3 复杂度分析

试除法的时间复杂度是 O(sqrt(n))。对于一个 1024 位的数，sqrt(n) 大约是 2^512——这个数字本身就有 154 位十进制数。即使每次除法只需要 1 纳秒，完成试除也需要大约 10^144 年。显然，这在实际应用中完全不可行。

1.4 优化：小素数筛

虽然试除法不能作为大数素性测试的主要手段，但它在工业级素数生成流水线中扮演了重要的预过滤角色。用前几千个小素数对候选数做试除，可以快速排除大量合数。统计表明，用前 2000 个素数做试除，可以排除大约 88% 的随机奇数。

二、Fermat 小定理与 Fermat 测试

2.1 定理陈述

Fermat 小定理是数论中最优美的结论之一：

若 p 是素数，且 gcd(a, p) = 1，则 a^(p-1) = 1 (mod p)。

逆否命题给出了一个素性测试方法：如果存在某个 a 使得 a^(n-1) != 1 (mod n)，则 n 一定是合数。这样的 a 被称为 n 的 Fermat 证人（witness）。

2.2 Fermat 测试的实现思路

/* Fermat 素性测试（伪代码风格） */
bool fermat_test(uint64_t n, int k)
{
    if (n < 2) return false;
    if (n < 4) return true;

    for (int i = 0; i < k; i++) {
        uint64_t a = 2 + rand() % (n - 3);  /* 随机选择 a in [2, n-2] */
        if (mod_pow(a, n - 1, n) != 1)
            return false;  /* n 一定是合数 */
    }
    return true;  /* n 可能是素数 */
}

2.3 Carmichael 数：Fermat 测试的致命缺陷

Fermat 测试看起来简单有效，但有一类特殊的合数能够欺骗它。Carmichael 数（也叫绝对伪素数）是满足以下条件的合数 n：对所有与 n 互素的 a，都有 a^(n-1) = 1 (mod n)。

最小的 Carmichael 数是 561 = 3 * 11 * 17。可以验证，对任意与 561 互素的 a，都有 a^560 = 1 (mod 561)。

Carmichael 数虽然稀少，但已被证明有无穷多个（1994 年由 Alford、Granville 和 Pomerance 证明）。前几个 Carmichael 数是：

561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, ...

Korselt 判别法给出了 Carmichael 数的充要条件：正合数 n 是 Carmichael 数当且仅当 n 是无平方因子的，且对 n 的每个素因子 p，都有 (p - 1) | (n - 1)。

2.4 为什么 Fermat 测试不够

Carmichael 数的存在意味着，无论进行多少轮 Fermat 测试，都无法将错误概率降到任意低——除非我们恰好选中了与 n 不互素的 a（此时直接得到了 n 的一个因子）。但对于大的 Carmichael 数，这种概率极低。

这就是为什么工业级实现不使用纯 Fermat 测试，而是使用它的加强版——Miller-Rabin 测试。

三、Miller-Rabin 测试：强伪素数

3.1 从 Fermat 到 Miller-Rabin

Miller-Rabin 测试的核心思想是利用二次探测定理对 Fermat 测试进行加强。

设 n 是奇素数，将 n - 1 写成 2^s * d 的形式，其中 d 是奇数。对于任意 1 <= a <= n-1，以下两个条件中至少有一个成立：

1. a^d = 1 (mod n)
2. 存在某个 r (0 <= r < s)，使得 a⁽²r * d) = -1 (mod n)

这来自于一个简单的代数事实：如果 x^2 = 1 (mod p) 且 p 是奇素数，则 x = 1 或 x = -1 (mod p)。换句话说，在模素数的运算中，1 的平方根只有 +1 和 -1。

3.2 算法描述

给定待测数 n 和基 a：

1. 将 n - 1 写成 2^s * d，d 为奇数
2. 计算 x = a^d mod n
3. 如果 x = 1 或 x = n - 1，则 n 对基 a 通过测试
4. 对 r = 1, 2, …, s-1：
   • x = x^2 mod n
   • 如果 x = n - 1，则 n 通过测试
   • 如果 x = 1，则 n 未通过测试（n 是合数）
5. 如果循环结束仍未返回，则 n 未通过测试

3.3 错误概率：4^(-k)

Miller-Rabin 测试最重要的性质是：对于任何奇合数 n，在 [1, n-1] 中随机选取基 a，n 通过测试的概率至多为 1/4。

这意味着，如果独立进行 k 轮测试，每轮随机选择一个基，那么一个合数通过所有 k 轮测试的概率至多为 (1/4)^k = 4^(-k)。

值得注意的是，4^(-k) 是最坏情况的上界。对于大多数合数，实际通过率远低于此。Damgard、Landrock 和 Pomerance 在 1993 年的论文中给出了更精确的分析，表明对于随机选择的大奇数，经过几轮 Miller-Rabin 测试后，误判概率远小于 4^(-k)。

3.4 为什么是 1/4 而不是 1/2

这个 1/4 的界来自于对强说谎者（strong liars）数量的精确分析。设 n 是奇合数，定义基 a 是 n 的强说谎者，如果 n 对基 a 通过 Miller-Rabin 测试。Monier 和 Rabin 独立证明了：强说谎者的数量至多为 phi(n)/4，其中 phi 是 Euler 函数。

直觉上，相比 Fermat 测试，Miller-Rabin 测试多了一层”二次探测”的检验，这使得说谎者的数量从 phi(n) 的一个较大比例降到了至多 1/4。

四、确定性 Miller-Rabin：小范围的精确测试

4.1 确定性见证集

虽然 Miller-Rabin 本质上是概率算法，但对于有限范围内的数，可以找到一组固定的基，使得测试变成确定性的。

已知的结果包括：

最后一行意味着，对于所有小于 3.3 * 10^24 的数，只需用 12 个固定的基进行 Miller-Rabin 测试，就能得到完全确定性的结果。

4.2 在广义 Riemann 猜想下的确定性

在广义 Riemann 猜想（GRH）成立的假设下，Miller 在 1976 年证明了：对于任意合数 n，在 2 到 2 * (ln n)^2 的范围内一定存在一个 Miller-Rabin 证人。这意味着只需测试 O((log n)^2) 个基，就能确定性地判定素性。

然而，GRH 至今未被证明，所以这个结果在理论上仍有条件。但在实践中，人们观察到所需的证人数量远小于这个理论上界。

4.3 竞赛编程中的应用

在算法竞赛中，经常需要对 64 位范围内的数做素性测试。以下 7 个基足以覆盖所有 64 位正整数：

{2, 3, 5, 7, 11, 13, 37}

更保守的选择是使用以下基，这已被验证对所有 64 位数正确：

{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}

五、完整的 Miller-Rabin C 实现

以下是一个完整的、可编译运行的 Miller-Rabin 素性测试实现，包含模幂运算和防溢出的模乘法。

/*
 * miller_rabin.c -- Miller-Rabin Primality Test
 *
 * 支持 64 位整数范围内的素性测试。
 * 使用 __uint128_t 防止模乘溢出。
 *
 * 编译: gcc -O2 -o miller_rabin miller_rabin.c
 * 用法: ./miller_rabin <number>
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <time.h>

/* ----------------------------------------------------------------
 * 模乘法：计算 (a * b) % m，使用 128 位中间结果防止溢出
 * ---------------------------------------------------------------- */
static inline uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t m)
{
    return (unsigned __int128)a * b % m;
}

/* ----------------------------------------------------------------
 * 模幂运算：计算 base^exp % mod，使用快速幂（二进制指数法）
 *
 * 时间复杂度: O(log exp) 次模乘
 * ---------------------------------------------------------------- */
static uint64_t pow_mod(uint64_t base, uint64_t exp, uint64_t mod)
{
    uint64_t result = 1;

    base %= mod;
    if (base == 0)
        return 0;

    while (exp > 0) {
        if (exp & 1)
            result = mul_mod(result, base, mod);
        exp >>= 1;
        base = mul_mod(base, base, mod);
    }
    return result;
}

/* ----------------------------------------------------------------
 * 单轮 Miller-Rabin 测试
 *
 * 给定 n-1 = 2^s * d（d 为奇数），以 a 为基进行测试。
 * 返回 true 表示 n 对基 a "可能是素数"（通过测试）。
 * 返回 false 表示 n 一定是合数。
 * ---------------------------------------------------------------- */
static bool miller_rabin_round(uint64_t n, uint64_t d, int s, uint64_t a)
{
    uint64_t x = pow_mod(a, d, n);

    if (x == 1 || x == n - 1)
        return true;

    for (int r = 1; r < s; r++) {
        x = mul_mod(x, x, n);
        if (x == n - 1)
            return true;
        if (x == 1)
            return false;
    }
    return false;
}

/* ----------------------------------------------------------------
 * Miller-Rabin 素性测试（确定性版本，适用于 64 位整数）
 *
 * 对于 n < 2^64，使用以下 12 个基可以保证正确性：
 *   {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}
 *
 * 这组基在 n < 3.3 * 10^24 范围内被证明是完备的，
 * 而 2^64 < 3.3 * 10^24，因此覆盖了整个 uint64_t 范围。
 * ---------------------------------------------------------------- */
bool is_prime(uint64_t n)
{
    /* 处理小数和偶数 */
    if (n < 2)
        return false;
    if (n < 4)
        return true;
    if (n % 2 == 0)
        return false;

    /* 小素数快速判定 */
    static const uint64_t small_primes[] = {
        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37
    };
    static const int n_small = sizeof(small_primes) / sizeof(small_primes[0]);

    for (int i = 0; i < n_small; i++) {
        if (n == small_primes[i])
            return true;
        if (n % small_primes[i] == 0)
            return false;
    }

    /* 将 n-1 分解为 2^s * d */
    uint64_t d = n - 1;
    int s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        s++;
    }

    /* 用每个基进行测试 */
    for (int i = 0; i < n_small; i++) {
        uint64_t a = small_primes[i];
        if (a >= n)
            continue;
        if (!miller_rabin_round(n, d, s, a))
            return false;
    }
    return true;
}

/* ----------------------------------------------------------------
 * 概率版 Miller-Rabin（用于演示随机基的版本）
 *
 * 使用 k 轮随机基进行测试。
 * 错误概率上界: 4^(-k)
 * ---------------------------------------------------------------- */
bool is_prime_probabilistic(uint64_t n, int k)
{
    if (n < 2)
        return false;
    if (n < 4)
        return true;
    if (n % 2 == 0)
        return false;

    uint64_t d = n - 1;
    int s = 0;
    while ((d & 1) == 0) {
        d >>= 1;
        s++;
    }

    for (int i = 0; i < k; i++) {
        /* 随机选择 a in [2, n-2] */
        uint64_t a = 2 + (uint64_t)rand() % (n - 3);
        if (!miller_rabin_round(n, d, s, a))
            return false;
    }
    return true;
}

/* ----------------------------------------------------------------
 * 素数生成：找到不小于 start 的下一个素数
 * ---------------------------------------------------------------- */
uint64_t next_prime(uint64_t start)
{
    if (start <= 2)
        return 2;
    if (start % 2 == 0)
        start++;

    while (!is_prime(start)) {
        start += 2;
    }
    return start;
}

/* ----------------------------------------------------------------
 * 基准测试：测量不同轮数下的测试时间
 * ---------------------------------------------------------------- */
static void benchmark(void)
{
    /* 一些已知的大素数用于测试 */
    static const uint64_t test_primes[] = {
        999999999999999989ULL,
        999999999999999877ULL,
        18446744073709551557ULL,  /* 最大的 64 位素数之一 */
    };
    static const int n_tests = sizeof(test_primes) / sizeof(test_primes[0]);

    printf("\n--- Benchmark: Miller-Rabin Performance ---\n");
    printf("%-22s  %8s  %8s  %8s\n",
           "Number", "k=1", "k=10", "k=64");
    printf("------------------------------------------------------------\n");

    int rounds[] = {1, 10, 64};

    for (int t = 0; t < n_tests; t++) {
        uint64_t n = test_primes[t];
        printf("%-22lu", (unsigned long)n);

        for (int ri = 0; ri < 3; ri++) {
            int k = rounds[ri];
            int iterations = 100000;
            clock_t start = clock();
            for (int i = 0; i < iterations; i++) {
                is_prime_probabilistic(n, k);
            }
            clock_t end = clock();
            double us = (double)(end - start) / CLOCKS_PER_SEC * 1e6
                        / iterations;
            printf("  %6.2f us", us);
        }
        printf("\n");
    }
}

/* ----------------------------------------------------------------
 * 主函数
 * ---------------------------------------------------------------- */
int main(int argc, char *argv[])
{
    srand((unsigned)time(NULL));

    if (argc < 2) {
        fprintf(stderr, "Usage: %s <number> [--bench]\n", argv[0]);
        return 1;
    }

    if (strcmp(argv[1], "--bench") == 0) {
        benchmark();
        return 0;
    }

    uint64_t n = strtoull(argv[1], NULL, 10);
    if (n == 0) {
        fprintf(stderr, "Invalid number: %s\n", argv[1]);
        return 1;
    }

    bool result = is_prime(n);
    printf("%lu is %s\n", (unsigned long)n,
           result ? "PRIME" : "COMPOSITE");

    if (result) {
        uint64_t np = next_prime(n + 1);
        printf("Next prime: %lu\n", (unsigned long)np);
    }

    return 0;
}

5.1 代码要点

模乘法。64 位乘法可能产生 128 位的中间结果。我们利用 GCC 扩展的 __uint128_t 类型来处理这个问题。在不支持 128 位整数的平台上，需要使用 Montgomery 乘法或 Barrett 约简等技术。

快速幂。模幂运算 a^e mod m 使用二进制指数法，将 O(e) 次乘法降到 O(log e) 次。每次乘法后立即取模，保证中间结果不超出 128 位。

确定性见证集。对于 64 位范围内的数，使用 {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37} 这 12 个基可以保证零误判。这避免了随机数生成器的质量问题。

六、Baillie-PSW 测试：无已知反例

6.1 测试描述

Baillie-PSW 测试是 Miller-Rabin 的进一步强化，由 Robert Baillie、Carl Pomerance、John Selfridge 和 Samuel Wagstaff Jr. 在 1980 年提出。它结合了两种不同类型的素性测试：

1. 以 2 为基的强伪素数测试（即一轮 Miller-Rabin，基为 2）
2. 强 Lucas 伪素数测试

这两种测试基于不同的数学结构，因此它们的”弱点”几乎不重叠。

6.2 核心思想

Miller-Rabin 基于费马小定理，本质上是在乘法群 (Z/nZ)* 中检验元素的阶。而 Lucas 测试基于 Lucas 序列的性质，本质上是在某种二次扩张中检验”特征根”的行为。

一个合数要同时欺骗这两种测试，需要在两种完全不同的代数结构中同时”表现得像素数”，这在实践中极难发生。

6.3 反例搜索现状

截至目前，没有人找到任何一个 Baillie-PSW 伪素数。搜索已经覆盖了 2^64 以下的所有数。Carl Pomerance 甚至悬赏 620 美元（后来提高到 30 美元——这是一个数学家式的幽默，他说最初的悬赏太高了）给找到反例的人。

严格来说，这不意味着反例不存在。理论分析表明 Baillie-PSW 伪素数可能存在，但密度极低——可能比 10^10000 还大。对于实际应用而言，这和”不存在”没有区别。

6.4 实际影响

许多数学软件和密码学库将 Baillie-PSW 作为标准素性测试：

• Mathematica 的 PrimeQ 函数
• Maple 的 isprime 函数
• GMP 库的 mpz_probab_prime_p（部分）
• Python 3.9+ 的 sympy.isprime

由于该测试只需要一轮 Miller-Rabin（基 2）加一轮 Lucas 测试，其速度接近于两轮 Miller-Rabin，但提供的可靠性远超任何有限轮数的 Miller-Rabin。

七、AKS 素性测试：理论上的多项式时间

7.1 历史意义

2002 年，印度计算机科学家 Manindra Agrawal、Neeraj Kayal 和 Nitin Saxena 发表了论文 “PRIMES is in P”，给出了第一个确定性的、无条件的多项式时间素性测试算法。这解决了一个长期悬而未决的理论问题，也使三位作者获得了 Godel 奖和 Fulkerson 奖。

7.2 算法思想

AKS 算法基于以下推广的费马小定理：

n 是素数当且仅当 (X + a)^n = X^n + a (mod n) 对所有 a 成立。

直接验证这个等式需要展开 n+1 个系数，计算量为 O(n)——对于大数来说太慢。AKS 的关键突破是将这个等式模一个精心选择的多项式 X^r - 1，从而将系数数量降到 O(r)，其中 r 是一个关于 log n 的多项式量级的数。

7.3 复杂度

AKS 算法的原始复杂度为 O((log n)^12)，后来被改进到 O((log n)^6)。在 Lenstra 和 Pomerance 的改进版本中，复杂度进一步降到约 O((log n)^6)。

7.4 实际表现

尽管 AKS 在理论上是多项式时间的，但其常数因子非常大。对于实际大小的数（比如 1024 位），AKS 比 Miller-Rabin 慢数个数量级。

粗略估计：

• Miller-Rabin（64 轮）对 1024 位数: 约 50 毫秒
• AKS 对 1024 位数: 可能需要数小时甚至数天

因此，AKS 在密码学实践中几乎不被使用。它的价值是理论上的：证明了素性测试问题属于复杂度类 P。

7.5 学术价值与工程价值的分离

AKS 算法是”理论上最优但实践中无用”的典型案例。这提醒我们：

• 渐近复杂度不等于实际性能
• 概率算法在实践中往往优于确定性算法
• 密码学工程需要在理论纯粹性和实际效率之间做取舍

八、Lucas-Lehmer 测试：Mersenne 素数的专用工具

8.1 Mersenne 素数

形如 M_p = 2^p - 1 的素数称为 Mersenne 素数。目前已知最大的素数几乎都是 Mersenne 素数，这要归功于 Lucas-Lehmer 测试的高效性。

截至 2024 年，已知最大的 Mersenne 素数是 M_136279841 = 2^136279841 - 1，它有超过 4100 万位十进制数。

8.2 Lucas-Lehmer 测试

对于 p > 2，M_p = 2^p - 1 是素数当且仅当 S_(p-2) = 0 (mod M_p)，其中 Lucas-Lehmer 序列定义为：

S_0 = 4
S_i = S_(i-1)^2 - 2 (mod M_p)    对 i >= 1

8.3 实现

/*
 * Lucas-Lehmer 测试（简化版，使用 GMP 大数库）
 * 实际实现需要链接 -lgmp
 */

#include <gmp.h>
#include <stdbool.h>

bool lucas_lehmer(unsigned long p)
{
    if (p == 2)
        return true;
    if (p % 2 == 0)
        return false;

    mpz_t mp, s, temp;
    mpz_init(mp);
    mpz_init(s);
    mpz_init(temp);

    /* M_p = 2^p - 1 */
    mpz_ui_pow_ui(mp, 2, p);
    mpz_sub_ui(mp, mp, 1);

    /* S_0 = 4 */
    mpz_set_ui(s, 4);

    /* 迭代 p-2 次 */
    for (unsigned long i = 0; i < p - 2; i++) {
        mpz_mul(s, s, s);       /* s = s^2 */
        mpz_sub_ui(s, s, 2);   /* s = s - 2 */
        mpz_mod(s, s, mp);     /* s = s mod M_p */
    }

    bool result = (mpz_cmp_ui(s, 0) == 0);

    mpz_clear(mp);
    mpz_clear(s);
    mpz_clear(temp);

    return result;
}

8.4 为什么 Lucas-Lehmer 如此高效

Lucas-Lehmer 测试只需要 p - 2 次模平方运算，每次运算涉及一个约 p 位的数。使用 FFT 乘法，每次模平方的时间为 O(p * log p * log log p)，总时间为 O(p^2 * log p * log log p)。

相比之下，通用的素性测试（如 Miller-Rabin）对 2^p - 1 的测试需要 O(p) 次模幂运算，每次模幂本身又需要 O(p) 次模乘，总复杂度更高。

Lucas-Lehmer 测试之所以高效，是因为它利用了 Mersenne 数的特殊代数结构。这也是为什么大素数搜索项目（如 GIMPS）专注于 Mersenne 素数。

九、工业级素数生成流水线

9.1 总体架构

工业级素数生成（如 RSA 密钥生成中的素数选取）通常采用以下流水线：

1. 随机候选数生成: 生成一个 n 位的随机奇数，设置最高位为 1（保证位长度）
2. 小素数试除: 用前几千个小素数做试除，快速排除大量合数
3. Miller-Rabin 测试: 对通过试除的候选数进行多轮 Miller-Rabin 测试
4. 输出素数: 通过所有测试的数作为素数输出

9.2 候选数生成策略

有两种常见的候选数生成策略：

随机跳跃法：每次生成一个全新的随机数。优点是各候选数之间独立，安全性分析简单。缺点是需要频繁调用随机数生成器。

增量搜索法：生成一个随机起点 n，然后依次测试 n, n+2, n+4, …（只测试奇数）。优点是后续候选数可以增量更新试除结果。缺点是候选数之间存在相关性，但实际分析表明这不会降低安全性。

OpenSSL 使用增量搜索法，这也是目前最主流的做法。

9.3 素数定理与期望搜索次数

根据素数定理，一个随机 n 位数是素数的概率约为 1/(n * ln 2)。因此，找到一个 n 位素数平均需要测试约 n * ln 2 个候选数。

9.4 小素数试除的详细实现

/*
 * 使用小素数筛进行预过滤。
 * 维护一个模残差表，支持增量更新。
 */

#define N_SMALL_PRIMES 2048

static const uint16_t small_primes[N_SMALL_PRIMES] = {
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
    /* ... 省略中间部分，实际使用时需要完整的素数表 ... */
    /* 第 2048 个素数是 17863 */
};

/*
 * 检查 n 是否能被任何小素数整除。
 * 返回 0 如果 n 通过试除（可能是素数），
 * 返回具体的小素因子如果 n 是合数。
 */
uint16_t trial_division(uint64_t n)
{
    for (int i = 0; i < N_SMALL_PRIMES; i++) {
        if (n % small_primes[i] == 0) {
            if (n == small_primes[i])
                return 0;  /* n 本身就是小素数 */
            return small_primes[i];
        }
    }
    return 0;
}

/*
 * 增量试除：已知 n 对各小素数的模残差，
 * 计算 n+delta 的残差只需一次加法和取模。
 */
typedef struct {
    uint16_t residues[N_SMALL_PRIMES];
} sieve_state_t;

void sieve_init(sieve_state_t *state, uint64_t n)
{
    for (int i = 0; i < N_SMALL_PRIMES; i++)
        state->residues[i] = n % small_primes[i];
}

bool sieve_check(const sieve_state_t *state)
{
    for (int i = 1; i < N_SMALL_PRIMES; i++) {
        if (state->residues[i] == 0)
            return false;  /* 能被 small_primes[i] 整除 */
    }
    return true;
}

void sieve_advance(sieve_state_t *state, uint16_t delta)
{
    for (int i = 0; i < N_SMALL_PRIMES; i++) {
        state->residues[i] = (state->residues[i] + delta) % small_primes[i];
    }
}

十、OpenSSL 的 BN_is_prime_ex 实现解析

10.1 函数签名

OpenSSL 中素性测试的核心函数是 BN_check_prime（较新版本）和旧版的 BN_is_prime_ex：

/* OpenSSL 3.x API */
int BN_check_prime(const BIGNUM *p, BN_CTX *ctx, BN_GENCB *cb);

/* OpenSSL 1.x 旧 API */
int BN_is_prime_ex(const BIGNUM *p, int nchecks,
                   BN_CTX *ctx, BN_GENCB *cb);

10.2 内部流程

以 OpenSSL 3.x 的 BN_check_prime 为例，其内部流程大致如下：

BN_check_prime(p)
  |
  +-> 步骤 1: 检查 p 是否为小素数或能被小素数整除
  |     使用 BN_mod_word 对前 2048 个素数做试除
  |
  +-> 步骤 2: 一轮 Miller-Rabin 测试，基为随机选择
  |     调用 BN_is_prime_fasttest_ex 的内部逻辑
  |
  +-> 步骤 3: 一轮 Lucas 测试
  |     实现 Baillie-PSW 测试的 Lucas 部分
  |
  +-> 步骤 4: 额外的 Miller-Rabin 轮次
  |     轮次数由安全级别决定，默认参考 FIPS 186-4
  |
  +-> 返回结果: 1 = 可能是素数, 0 = 合数, -1 = 错误

10.3 轮次选择策略

OpenSSL 根据待测数的位长度自动选择 Miller-Rabin 的轮次数。这遵循 FIPS 186-4 的建议：

注意，随着位长度增加，所需的 Miller-Rabin 轮次反而减少。这是因为对于更大的数，随机合数通过 Miller-Rabin 的概率更低（Damgard-Landrock-Pomerance 的结果）。

10.4 回调机制

OpenSSL 的素数生成支持回调机制，允许应用程序监控生成进度：

/*
 * 回调函数在以下时刻被调用:
 * - 生成新的候选数时 (event = 0, ret = '.')
 * - 候选数通过试除时 (event = 1, ret = '+')
 * - 候选数通过 Miller-Rabin 时 (event = 2, ret = '*')
 *
 * 这就是为什么 openssl genrsa 命令会输出 '.', '+', '*' 等符号。
 */
int progress_callback(int event, int n, BN_GENCB *cb)
{
    char c = '.';
    if (event == 1) c = '+';
    if (event == 2) c = '*';
    putc(c, stderr);
    fflush(stderr);
    return 1;  /* 返回 0 会取消生成 */
}

如果你用过 openssl genrsa 命令，可能见过类似这样的输出：

Generating RSA private key, 2048 bit long modulus
...........++++
.......++++

其中每个 . 表示一个被试除排除的候选数，+ 表示通过了一轮 Miller-Rabin，++++ 表示通过了多轮测试。

十一、实际应用场景

11.1 RSA 密钥生成

RSA 密钥生成的核心步骤：

1. 生成两个大素数 p 和 q（各为密钥长度的一半，如 RSA-2048 需要两个 1024 位素数）
2. 计算 n = p * q
3. 计算 phi(n) = (p-1)(q-1)
4. 选择公钥指数 e（通常为 65537）
5. 计算私钥指数 d = e^(-1) mod phi(n)

在步骤 1 中，对素数 p 和 q 有额外的要求：

安全要求:
- |p - q| 应该足够大（防止 Fermat 因式分解）
- p - 1 和 q - 1 应该有大的素因子（防止 Pollard p-1 攻击）
- p 和 q 不应太接近 sqrt(n)
- FIPS 186-4 要求 |p - q| > 2^(nlen/2 - 100)

11.2 安全素数与 Diffie-Hellman

Diffie-Hellman 密钥交换通常需要安全素数（safe primes）。一个素数 p 是安全素数，当且仅当 (p-1)/2 也是素数。(p-1)/2 称为 Sophie Germain 素数。

安全素数的生成比普通素数更耗时，因为需要同时验证 p 和 (p-1)/2 的素性。

/*
 * 安全素数生成策略:
 * 1. 生成候选 Sophie Germain 素数 q
 * 2. 计算 p = 2q + 1
 * 3. 分别验证 q 和 p 的素性
 */
uint64_t generate_safe_prime(int bits)
{
    while (1) {
        uint64_t q = random_odd(bits - 1);  /* 生成 (bits-1) 位奇数 */

        if (!is_prime(q))
            continue;

        uint64_t p = 2 * q + 1;

        if (is_prime(p))
            return p;
    }
}

安全素数的密度大约是普通素数密度的 1/(ln n) 倍，因此生成时间大约是普通素数的 ln(n) 倍。对于 2048 位安全素数，平均需要测试约 (710)^2 / 2 ~ 250000 个候选数。

11.3 素数在椭圆曲线密码学中的角色

椭圆曲线密码学（ECC）中，素数主要出现在以下场景：

• 有限域 GF(p) 的模数 p 必须是素数
• 椭圆曲线群的阶通常要求是素数或接近素数
• 基点的阶必须是大素数

ECC 使用的素数通常有特殊形式（如 NIST P-256 使用的 Mersenne-like 素数），以加速模运算。

11.4 基准测试数据

以下是在不同硬件上，使用 OpenSSL 生成不同大小 RSA 密钥的基准时间：

时间随密钥长度超线性增长，主要原因是：

1. 素数密度降低：期望候选数数量线性增长
2. 模幂运算变慢：时间与位长度的立方成正比
3. 两者叠加导致总时间约与位长度的四次方成正比

十二、工程陷阱与个人观点

12.1 工程陷阱总览

以下总结了素性测试与素数生成中常见的工程陷阱：

12.2 侧信道攻击

素数生成过程中的侧信道泄露是一个严重但经常被忽视的问题。

时间侧信道。Miller-Rabin 测试的执行时间取决于候选数的值和基的选择。如果攻击者能够精确测量每次测试的耗时，可能推断出关于候选数的信息。

缓存侧信道。模幂运算中的平方-乘法模式取决于指数的二进制表示。如果攻击者能监控 CPU 缓存访问模式（如通过 Flush+Reload 攻击），可能恢复指数信息。

防御措施包括：

• 使用恒定时间的模幂运算（如 Montgomery 阶梯法）
• 在关键操作前后刷新缓存
• 使用操作系统提供的安全内存分配（如 mlock）

12.3 关于测试轮次的选择

工程实践中，Miller-Rabin 轮次的选择需要在安全性和性能之间取舍。我的建议是：

对于 RSA 密钥生成：遵循 FIPS 186-4 的建议就足够了。具体来说，对于 2048 位 RSA，4 轮 Miller-Rabin 加一轮 Lucas 测试（即 Baillie-PSW）提供的安全级别已远超 128 位。没有必要使用 64 轮或更多。

对于安全素数生成：由于安全素数的生成本身就很慢，额外几轮 Miller-Rabin 的开销相对较小。建议对 q 和 p 分别做完整的 Baillie-PSW 测试。

对于竞赛编程：如果数在 64 位范围内，使用确定性见证集即可，无需随机基。如果数更大（如 Python 的任意精度整数），5-10 轮随机基的 Miller-Rabin 通常足够。

12.4 概率算法的哲学思考

Miller-Rabin 测试是概率算法在工程实践中大获成功的典型案例。它带来了一个有趣的哲学问题：我们能信任概率算法吗？

答案是肯定的，而且理由很充分。64 轮 Miller-Rabin 的误判概率为 4^(-64) = 2^(-128)。作为对比：

• 同一个人连续被闪电击中 10 次的概率：约 10^(-70)
• 宇宙中所有原子同时发生量子隧穿的概率：远大于 2^(-128)
• 计算机硬件发生未检测到的位翻转的概率：约 10^(-15) 每小时

换句话说，你的计算机因为宇宙射线导致的位翻转而给出错误结果的概率，远高于 Miller-Rabin 将合数误判为素数的概率。确定性算法也会因硬件错误而给出错误结果。从这个角度看，Miller-Rabin 测试比任何确定性算法都”更可靠”——前提是你信任概率论。

12.5 关于 AKS 与实用主义

AKS 算法的发现是数学上的重大突破，但它在工程界几乎无人使用。这并不意味着 AKS 没有价值——它的价值在于：

1. 解决了 “PRIMES in P?” 这个长期开放问题
2. 带动了相关数论和算法研究的发展
3. 启发了对其他问题的去随机化（derandomization）研究

但在工程层面，我们应该务实地选择工具。Miller-Rabin 的 2^(-128) 误判率在工程上等价于零，而它的速度比 AKS 快几个数量级。正如 Donald Knuth 所说：

“过早优化是万恶之源，但过早追求理论完美同样是一种浪费。”

（好吧，后半句是我编的。但精神是一致的。）

12.6 未来展望

随着量子计算的发展，RSA 和基于大整数分解的密码方案面临威胁。Shor 算法可以在量子多项式时间内分解大整数，使 RSA 不再安全。然而，素性测试本身不会因此失去意义：

• 后量子密码学（如格密码、哈希签名）中仍然需要素数
• 椭圆曲线的构造需要素数阶的有限域
• 纠错码和其他代数结构中大量使用素数

素性测试作为计算数论的基础工具，其重要性将持续存在。

附录A：Miller-Rabin 测试的正确性证明概要

定理：若 n 是奇素数，则对任意 1 <= a <= n-1，n 对基 a 通过 Miller-Rabin 测试。

证明概要：

设 n - 1 = 2^s * d，其中 d 为奇数。

考虑序列：a^d, a^(2d), a^(4d), …, a⁽²s * d) = a^(n-1)

由 Fermat 小定理，a^(n-1) = 1 (mod n)。

如果 a^d = 1 (mod n)，则条件 1 成立，证毕。

否则，序列中存在第一个位置 r 使得 a⁽²r * d) = 1 (mod n) 但 a⁽²(r-1) * d) != 1 (mod n)。

设 x = a⁽²(r-1) * d)，则 x^2 = 1 (mod n) 但 x != 1 (mod n)。

因为 n 是素数，Z/nZ 是域，x^2 - 1 = (x-1)(x+1) = 0 (mod n) 意味着 x = 1 或 x = -1 (mod n)。

由于 x != 1，必有 x = n - 1 (mod n)，即条件 2 成立。

定理：对于任意奇合数 n，在 [1, n-1] 中至多有 (n-1)/4 个强说谎者。

这个定理的证明较为复杂，需要分别讨论 n 是素数幂、两个素数之积、以及一般情况。核心思想是利用中国剩余定理将 (Z/nZ)* 分解为分量，然后分析在每个分量中满足条件的元素数量。完整证明可参考 Rabin 的原始论文或 Monier 的独立证明。

附录B：相关算法复杂度对比

注：Baillie-PSW 标注”概率性“是因为理论上存在反例的可能，但至今未找到。ECPP 标注”确定性”是因为其运行时间分析是启发式的。

附录C：参考文献

1. Miller, G. L. “Riemann’s Hypothesis and Tests for Primality.” Journal of Computer and System Sciences, 1976.
2. Rabin, M. O. “Probabilistic Algorithm for Testing Primality.” Journal of Number Theory, 1980.
3. Baillie, R. and Wagstaff, S. S. “Lucas Pseudoprimes.” Mathematics of Computation, 1980.
4. Agrawal, M., Kayal, N., and Saxena, N. “PRIMES is in P.” Annals of Mathematics, 2004.
5. Damgard, I., Landrock, P., and Pomerance, C. “Average Case Error Estimates for the Strong Probable Prime Test.” Mathematics of Computation, 1993.
6. Alford, W. R., Granville, A., and Pomerance, C. “There Are Infinitely Many Carmichael Numbers.” Annals of Mathematics, 1994.
7. FIPS 186-4. “Digital Signature Standard (DSS).” NIST, 2013.
8. OpenSSL Source Code. crypto/bn/bn_prime.c.

上一篇: 快速傅里叶变换 下一篇: 扩展欧几里得与模逆元 相关阅读: - 椭圆曲线算术 - 有限域算术