有人可以帮我理解这个 Eratosthenes 脚本的筛子吗？最后几行让我难过

Question

好的，所以我明白了埃拉托色尼筛背后的原理。我尝试自己编写一个，但基本上失败了（我编写了一个功能素数生成器，它只是效率不高或筛子）。我认为我在理解数学方面遇到了更多问题，但是编程也让我感到困惑。

import numpy as np

def primesfrom2to(n):
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool) 
    # print(sieve) for n = 10 returns [True, True, True] 

好吧，算了，我有点糊涂了。它正在生成一个真值列表，当它们被确定为复合时，这些值将逐渐被标记为假。我认为这是说小于 n 的值中不超过 1/3 将是素数，但是添加模数运算可以实现什么？

    sieve[0] = False

将 1 标记为假？

    for i in range(int(n**0.5)//3+1):
        # print(i) for n = 10 returns 0 and 1 on separate lines

这是设置要检查的范围。没有数字的乘积大于其平方根。除以3+1 是怎么回事？

        if sieve[i]:
            k=3*i+1|1
            #print(k) for n = 10 returns '5' doesn't change till it adds 7 at n = 50

好的，如果sieve[i] 为真（所以素数/尚未标记为复合数？）那么3*i + 1 or 1？这究竟是如何工作的？似乎正在生成素数，稍后将相乘以删除所有后续产品，但我不知道如何。

            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False

好的，所以这两个都取素数k 并标记它们的所有进一步倍数？如果n=5 不是sieve[8.33::10] = False？另一个像sieve[41.3::10]？我只是不明白。

print(np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)])

好的，最后实际生成素数列表。再一次，将它乘以三是怎么回事？显然，这段代码中有一些关于 3 的基本内容，我只是不明白。我也无法再次掌握|1 的概念。

哦，只是为了好玩，这是我有效且低效的尝试。

import numpy
def sieve(num):

    master = numpy.array([i for i in range(2, num+1)])
    run = []
    y=2


    while y <= ((num+1)**(1/2)): 
        thing = [x for x in range(y, num+1, y) if x > 5 or x == 4]
        run = run + thing
        y = y+1 
    alist = numpy.in1d(master, run, invert = True)
    blist = (master[alist])
    print(blist)

计算素数到 500,000 需要 57 秒。我正在做 2,000,000 以内的欧拉素数和问题。

Answer 1

这是一个简单的 numpy 筛算法：

import numpy as np
def sieve(Nmax):
    """Calculate prime numbers between 0 and Nmax-1."""
    is_prime = np.ones(Nmax, dtype=bool)
    Pmax = int(np.sqrt(Nmax)+1)
    is_prime[0] = is_prime[1] = False
    p = 2
    dp = 1
    while p <= Pmax:
        if is_prime[p]:
            is_prime[p*2::p] = False
        p += dp
        dp = 2
    primes = np.argwhere(is_prime)[:,0]
    print(primes)

sieve(100)    

如果我删除打印语句，它会在我普通的笔记本电脑上在 2.5 秒内运行 Nmax=10^8。

这个算法效率有点低，因为它需要存储偶数和 3 的倍数的素数。你可以通过过滤掉它们来节省存储空间，这样你就可以跟踪 n* 的素数6+1 和 n*6+5，对于任何整数 n>=1。这为您节省了三分之二的存储空间，但代价是更多的簿记。您似乎试图从困难的版本开始。此外，如果涉及 Python 解释器评估 if 语句并进行列表的内存管理，所有的簿记往往会变得昂贵。 Python/numpy 的数组切片是一种更自然的方法，它可能对您的目的来说足够快。

关于您的问题：

• (n%6==2) 是布尔值，将被解释为 0 或 1，并且会再添加一个元素以防止出现“逐一”错误。
• int(n**0.5)//3+1 应读作int(int(n**0.5)/3) + 1。除法先于加法。除以三是为了只为 6n+1 和 6n+5 值分配空间。
• k=3*i+1|1 表示k=3*i+1，如果是偶数就加一个（按位或）。数组元素“i”对应于潜在的素数(3*i+1)|1。所以如果数组索引是[0, 1, 2, 3, 4, 5, ...]，它们代表数字[1, 5, 7, 11, 13, 17, ...]。
• 对于代表 6n+1 和 6n+5 的元素，将筛分元素设置为 False 是分别完成的，这就是为什么有两个赋值。
• 如果你在 Python 2.7 中运行它，整数除法总是被截断，所以 9/2 将得到 4。在 Python 3 中，最好使用 // 运算符进行整数除法，即 (k*k)/3应该是(k*k)//3。

编辑：您可能希望将您询问的算法归因于：Fastest way to list all primes below N。