本文详解如何修复基于 willans 公式的 python 素数生成器中因阶乘爆炸导致的 `overflowerror`，通过数学简化（利用余弦函数周期性）与高精度数值计算（`decimal` 模块）双重策略，实现稳定计算第 8 个及以上素数。
Willans 公式是一类利用三角函数和阶乘构造的素数判定/生成公式，其理论形式优美，但直接数值实现极易失败——核心症结在于：当 j 增大时，(factorial(j - 1) + 1) / j 迅速增长为超大有理数，而 math.cos() 要求浮点输入，factorial(7!) = 5040 尚可，但 factorial(10!) = 3628800 已逼近 float 精度极限；至 j=15 时，factorial(14) 超过 10¹⁰，强制转 float 必然触发 OverflowError。
关键洞察在于：cos(x) 是周期为 2π 的函数，即 cos(x) = cos(x mod 2π)。因此无需计算巨大实数 x = π × (factorial(j−1)+1)/j，而应先将其对 2π 取模，再代入 cos。但由于 factorial(j−1) 极大，直接计算 x % (2π) 仍会因中间值过大而失败。更稳健的做法是：利用模运算性质，将整个分数在模 2 意义下化简角度系数（因为 cos(π × r) = cos(π × (r mod 2))），即只保留 (factorial(j−1)+1)/j 的小数部分对 2 的余数。
然而，对任意大整数 a/b 高效计算 (a/b) mod 2 仍需避免除法溢出。此时 decimal 模块成为必要工具——它支持用户自定义精度的大数小数运算。以下是修复后的完整实现：
from decimal import Decimal, getcontext
import math
# 设置足够精度（例如 50 位小数，应对 n=10+）
getcontext().prec = 50
def nth_prime(n):
if not (isinstance(n, int) and n > 0):
raise ValueError("n must be a positive integer")
# Willans 公式核心：π(j) = Σ_{i=1}^j [cos²(π·( (i−1)!+1 )/i)]
# 其中 [·] 为 Iverson bracket（真为1，假为0），等价于 floor(cos²(...))
def is_prime_indicator(j):
if j == 1:
return 0 # 1 不是素数
# 计算 (factorial(j-1) + 1) / j，用 Decimal 避免 float 溢出
num = math.factorial(j - 1) + 1
denom = j
# 高精度除法
ratio = Decimal(num) / De
cimal(denom)
# 利用 cos(π·x) = cos(π·(x mod 2))，取小数部分对 2 的余数
x_mod2 = ratio % Decimal(2)
# 计算 cos(π * x_mod2)，转换为 float（此时 x_mod2 ∈ [0,2)，安全）
angle = float(x_mod2 * Decimal(math.pi))
cos_val = math.cos(angle)
return int(math.floor(cos_val ** 2 + 1e-15)) # 加小量防浮点截断误差
def prime_counting(i):
return sum(is_prime_indicator(j) for j in range(1, i + 1))
# Willans 的 nth prime 公式：p_n = 1 + Σ_{i=1}^{2^n} ⌊(n / π(i))^(1/n)⌋
end_sum = 0
upper_bound = 2 ** n
for i in range(1, upper_bound + 1):
pi_i = prime_counting(i)
if pi_i == 0:
continue # 避免除零
# 计算 (n / pi_i)^(1/n)，用 Decimal 保障精度
base = Decimal(n) / Decimal(pi_i)
root = base ** (Decimal(1) / Decimal(n))
end_sum += int(math.floor(float(root) + 1e-12))
return end_sum + 1
# 测试
print(nth_prime(1)) # 2
print(nth_prime(5)) # 11
print(nth_prime(8)) # 19 ✅ 不再溢出
cimal(denom)
# 利用 cos(π·x) = cos(π·(x mod 2))，取小数部分对 2 的余数
x_mod2 = ratio % Decimal(2)
# 计算 cos(π * x_mod2)，转换为 float（此时 x_mod2 ∈ [0,2)，安全）
angle = float(x_mod2 * Decimal(math.pi))
cos_val = math.cos(angle)
return int(math.floor(cos_val ** 2 + 1e-15)) # 加小量防浮点截断误差
def prime_counting(i):
return sum(is_prime_indicator(j) for j in range(1, i + 1))
# Willans 的 nth prime 公式：p_n = 1 + Σ_{i=1}^{2^n} ⌊(n / π(i))^(1/n)⌋
end_sum = 0
upper_bound = 2 ** n
for i in range(1, upper_bound + 1):
pi_i = prime_counting(i)
if pi_i == 0:
continue # 避免除零
# 计算 (n / pi_i)^(1/n)，用 Decimal 保障精度
base = Decimal(n) / Decimal(pi_i)
root = base ** (Decimal(1) / Decimal(n))
end_sum += int(math.floor(float(root) + 1e-12))
return end_sum + 1
# 测试
print(nth_prime(1)) # 2
print(nth_prime(5)) # 11
print(nth_prime(8)) # 19 ✅ 不再溢出注意事项与优化建议：
- 精度权衡：getcontext().prec 设置过高会降低性能，建议根据 n 动态调整（如 n≤10 时设 prec=30，n≤15 时设 prec=60）；
- 效率警示：Willans 公式时间复杂度为 O(2ⁿ × n!)，仅适用于教学或极小 n（n≤12），生产环境请使用埃氏筛或分段筛；
- 数值稳定性：cos²(x) 在 x 接近整数时接近 1 或 0，浮点误差可能误判，代码中添加了 1e-15 补偿；
- 替代方案：若仅需验证公式逻辑，可用 sympy 符号计算 cos(pi * Rational(a,b))，完全规避浮点误差。
综上，修复本质是将不可行的“大数→浮点→三角函数”链，重构为“大数→高精度有理/小数→模约简→安全浮点三角”流程。这不仅解决了溢出，更体现了数值计算中“数学简化优先于蛮力精度”的工程哲学。
