优化精确素数定理（EPNT）的计算速度

先看个实际例子：给你前499个素数的序列：

2,3,5,7,...,3541,3547,3557,3559 

要预测下一个素数的话，答案是3571——也就是第500个素数。

素数定理（PNT） #

素数定理（PNT）能快速估算第n个素数，对应的近似公式是：

用它计算p₅₀₀ ≈ 3107，只需要微秒级的时间，速度快得离谱！

精确素数定理（EPNT） #

我做了个实验性的精确素数定理（EPNT），它能算出完全精确的第n个素数，公式是：

但尴尬的是，用它算p₅₀₀ = 3571居然要花25分钟！

问题来了 #

目前这个EPNT已经能准确算出前500个素数，但要验证更大的素数时，速度慢得让人崩溃！有没有优化技巧能提升它的计算速度？我自己先想到了几个方向：

• 换用Python以外的编程语言
• 引入多线程并行计算
• 替换为更快的高精度数学库
• 运行时动态调整小数精度参数mp.dps
• 借助WolframAlpha这类专业数学计算引擎

下面是我目前在用的Python实现代码：

import time
from mpmath import ceil, ln, mp, mpf, exp, fsum, power, zeta
from sympy import symbols, Eq, pprint, prime

N=500   # <--- 计算第N个素数
    
mp.dps = 20000
primes = []

def vengy_prime(k):
    # 确定性计算第k个素数
    s = ceil(k * ln(k * ln(k)))
    
    # 确定动态的Rosser(1941)上界
    N = int(ceil(k * (ln(k) + ln(ln(k)))))
    
    # 计算到N的有限求和
    print(f"Computing {N} zeta terms ...")   
    start_time = time.time()
    sum_N = fsum([1 / power(mpf(n), s) for n in range(1, N)])
    end_time = time.time()      
    print(f"Time taken: {end_time - start_time:.6f} seconds")
    
    # 计算涉及前k-1个素数的乘积项
    print(f"Computing product of {k-1} previous primes ...")  
    start_time = time.time()
    prod = exp(fsum([ln(1 - power(p, -s)) for p in primes[:k-1]]))
    end_time = time.time()      
    print(f"Time taken: {end_time - start_time:.6f} seconds")
    
    # 计算下一个素数p_k
    p_k=ceil((1 - 1 / (sum_N * prod)) ** (-1 / s))
    return p_k


# 生成已知的前k-1个素数
print("\nListing", N-1, "known primes:")
for k in range(1, N):
    p = prime(k)
    primes.append(p)        
print(primes)        

primes.append(vengy_prime(N))
pprint(Eq(symbols(f'p_{N}'), int(primes[-1])))

最新进展：优化后速度直接起飞！ #

多亏了Jérôme Richard帮忙优化的代码，现在算第500个素数只需要10秒！
优化后的运行日志如下：

Computing 4021 zeta terms ...
Time taken: 7.968423 seconds
Computing product of 499 previous primes ...
Time taken: 1.960771 seconds
p₅₀₀ = 3571

对比我原来的代码，算同样的内容居然要1486秒：

Computing 4021 zeta terms ...
Time taken: 1173.899538 seconds
Computing product of 499 previous primes ...
Time taken: 313.833039 seconds
p₅₀₀ = 3571

更惊喜的是，优化后的代码计算第4000个素数（参数设置N = 4000, precision = 700000），也只花了45分钟！

备注：内容来源于stack exchange，提问作者vengy