素数定理：用Python实现Li(x)积分以估算大区间素数数量

1. Li(x)积分用于素数数量估算是否正确？ #

• 完全正确。素数定理明确指出，当x趋近于无穷大时，素数计数函数π(x)（即≤x的素数个数）与对数积分Li(x)的比值趋近于1。相比x/ln(x)，Li(x)的近似精度更高，尤其是针对10^14这类超大数值，其误差远小于前者，是更优的素数数量近似方法。

2. Python中如何实现Li(x)积分？ #

对数积分Li(x)的常用定义是取积分下限为2（因ln(1)无意义，且小于2的素数仅1个，不影响大x的近似效果），即Li(x) = ∫₂ˣ dt/ln(t)。以下是两种可行的实现方式：

方法一：使用SciPy的高精度积分函数（推荐） #

SciPy的quad函数能高效处理反常积分，精度足够应对超大x的计算需求：

from scipy.integrate import quad
import numpy as np

def li(x):
    if x <= 2:
        return 0.0  # 若需严格对应π(2)=1，可单独判断返回1
    result, error = quad(lambda t: 1 / np.log(t), 2, x)
    return result

# 示例：计算10^14的对数积分
x = 10**14
approx = li(x)
print(f"Li({x}) ≈ {approx:.2f}")

对于超大x，quad的计算速度快且误差极小；若需适配小x的素数计数，可在函数内添加判断逻辑。

方法二：手动实现数值积分（无依赖场景） #

若无法使用SciPy，可采用辛普森法则结合对数采样实现，避免大x时区间分布不均：

import numpy as np

def li_manual(x, num_points=100000):
    if x <= 2:
        return 0.0
    # 对数采样让小t区间的点更密集，保证积分精度
    t = np.logspace(np.log10(2), np.log10(x), num_points)
    dt = np.diff(t)
    y = 1 / np.log(t)
    # 辛普森法则计算积分
    integral = np.sum((y[:-1] + y[1:]) * dt / 2)
    return integral

# 示例
x = 10**14
approx = li_manual(x)
print(f"Li({x}) ≈ {approx:.2f}")

可通过调整num_points参数平衡精度与计算速度，点数越多精度越高，但耗时也会增加。

内容的提问来源于stack exchange，提问作者Arthur Brenno