π(x)素数计数代码问题排查：N≤1e5正常、1000000结果错误

以下代码用于计算小于等于N的素数个数，在N≤100000的范围内运行结果完全准确：

• 输入输出对照表

|    输入    |   代码输出      | 
|-------------|---------------|
| 10          |     4✔       |
| 100         |     25✔      |
| 1000        |     168✔     |
| 10000       |     1229✔    |
| 100000      |     9592✔    |
| 1000000     |     78521✘   | 

而π(1000000)的正确值为78498，可见代码在N超过100000时计算结果存在误差。

原实现代码：

import time
def pi(x):
    nums = set(range(3,x+1,2))
    nums.add(2)
    #print(nums)
    prm_lst = set([])
    while nums:
        p = nums.pop()
        prm_lst.add(p)
        nums.difference_update(set(range(p, x+1, p)))
        #print(prm_lst)
    return prm_lst


if __name__ == "__main__":
    N = int(input())
    start = time.time()
    print(len(pi(N)))
    end= time.time()
    print(end-start)

问题分析 #

核心错误原因 #

这段代码本意是实现埃拉托斯特尼筛法，但存在致命的隐式依赖问题：Python的set是无序容器，pop()方法没有固定的弹出顺序，仅在待处理元素数值小、数量少的时候，因为底层哈希表无冲突，刚好按从小到大的顺序弹出元素，此时每次弹出的都是未被筛除的最小素数，逻辑和标准埃氏筛一致，所以小范围N计算正确。
当N增长到1000000量级时，set存储的整数数量变多，哈希冲突概率上升，pop()的弹出顺序不再是升序，会优先弹出哈希表中靠前的合数，代码会误将这个合数判定为素数加入结果集合，再筛除它的倍数，最终导致素数统计结果偏大。

优化方案 #

直接采用标准埃氏筛的稳定实现，用布尔数组代替set存储素数标记，完全规避容器顺序依赖问题，同时运行效率和内存利用率都远高于原实现：

import time
import math

def pi(x):
    if x < 2:
        return 0
    is_prime = [True] * (x + 1)
    is_prime[0] = is_prime[1] = False
    for i in range(2, int(math.isqrt(x)) + 1):
        if is_prime[i]:
            # 从i*i开始标记倍数，前面的已经被更小的素数标记过，减少重复操作
            is_prime[i*i : x+1 : i] = [False] * len(is_prime[i*i : x+1 : i])
    return sum(is_prime)

if __name__ == "__main__":
    N = int(input())
    start = time.time()
    print(pi(N))
    end = time.time()
    print(f"耗时: {end - start}s")

如果需要计算更大数值的π(x)（比如1e10以上），全量筛法会占用大量内存，可以改用专门用于素数计数的Meissel-Lehmer算法，仅需计算部分区间的素数就能得到准确结果，内存开销和运行速度都有量级提升。

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