Willans 公式实现素数生成时的数值溢出解决方案_技术教程

本文介绍如何修复基于 willans 公式（利用三角函数与阶乘判断素数）的 python 素数生成代码中因阶乘过大导致 `overflowerror` 的问题，核心是避免将超大整数直接转为浮点数计算余弦值。

Willans 公式是一类基于初等函数（如三角函数、取整函数）构造的“显式”素数公式，其理论形式优美，但在实际编程中极不实用——尤其当涉及 factorial(j-1) 时，j 仅需超过 20 就会使阶乘值远超 IEEE 754 双精度浮点数的表示范围（约 10^308），而 cos() 函数要求输入为 float，强制类型转换即触发 OverflowError。

你遇到的错误：

sum += floor(pow(cos(pi * (factorial(j - 1) + 1) / j), 2))

根本原因并非 cos 本身，而是 (factorial(j - 1) + 1) / j 这一除法在 Python 中默认产生 float（即使分子是超大整数），而 factorial(7) 已达 5040，factorial(15) 超过 1.3×10¹²，到 j=25 时 factorial(24) ≈ 6.2×10²³，已超出 float 精度并逼近溢出阈值——这正是 nth_prime(8) 在内部循环中 j 增大后崩溃的根源。

⚠️ 关键误区：decimal 模块无法解决此问题，因为 math.cos() 不接受 Decimal 类型，且三角函数库依赖底层 C double 实现，无法绕过浮点限制。

✅ 正确解法：利用余弦函数的周期性——cos(x) = cos(x mod 2π)。我们只需将大数参数对 2π 取模，再传入 cos。但由于 π 是无理数，直接模运算仍需高精度浮点。更稳健、可落地的做法是：用数论性质替代三角计算。

Willans 公式中关键判别项：

\left\lfloor \cos^2\left(\pi \cdot \frac{(j-1)! + 1}{j}\right) \right\rfloor
= 
\begin{cases}
1 & \text{if } j \text{ is prime (by Wilson's Theorem)} \\
0 & \text{if } j > 4 \text{ is composite}
\end{cases}

而 Wilson 定理指出：j 是素数 ⇔ (j-1)! ≡ -1 (mod j) ⇔ (j-1)! + 1 ≡ 0 (mod j) ⇔ j 整除 (j-1)! + 1。

因此，原式中的 cos²(...) 实质是在判断 j | ((j-1)! + 1)。我们完全可用纯整数模运算替代：

def is_prime_by_wilson(j):
    if j < 2:
        return False
    if j == 2:
        return True
    if j % 2 == 0:
        return False
    # 计算 (j-1)! mod j —— 注意：若 j 是合数且非平方因子，(j-1)! ≡ 0 (mod j)
    # 但为严谨，直接计算 (j-1)! + 1 是否能被 j 整除（仅对小 j 可行）
    try:
        fact_mod_j = 1
        for i in range(2, j):
            fact_mod_j = (fact_mod_j * i) % j
        return (fact_mod_j + 1) % j == 0
    except OverflowError:
        return False  # 防御性返回（实际 j 较小时不会触发）

# 但注意：Wilson 判定本身时间复杂度 O(j)，不适用于大 j；此处仅用于教学替换思路

然而，更现实的工程建议是：放弃 Willans 公式用于实际素数生成。它的时间复杂度为 O(2ⁿ × n!)，空间与数值稳定性均不可控。推荐改用高效可靠的方法，例如：

✅ 埃氏筛法（Eratosthenes）：生成前 N 个素数（预估上界用 n(log n + log log n)）
✅ 分段筛 + 计数优化：适合求第 n 个素数
✅ Miller-Rabin + 试除回溯：对大 n 高效验证

示例：用筛法稳健求第 8 个素数（即 19）：

def nth_prime_sieve(n):
    if n < 1:
        raise ValueError("n must be positive")
    # 保守估计第 n 个素数的上界（n>=6 时成立）
    limit = max(12, int(n * (log(n) + log(log(n))))) + 10
    sieve = [True] * (limit + 1)
    sieve[0] = sieve[1] = False
    for i in range(2, int(limit**0.5) + 1):
        if sieve[i]:
            sieve[i*i : limit+1 : i] = [False] * len(sieve[i*i : limit+1 : i])
    primes = [i for i, is_p in enumerate(sieve) if is_p]
    if len(primes) < n:
        return nth_prime_sieve(n + 1)  # 动态扩容（生产环境应更优估计）
    return primes[n-1]

from math import log
print(nth_prime_sieve(8))  # 输出: 19

? 总结：