有人可以帮助我理解这个埃拉托色尼筛法脚本吗？最后几行让我难过

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]

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

    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)

计算出50万个素数用了57s。我在做欧拉素数和高达 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 的倍数的质数。您可以通过过滤掉这些来节省存储空间 space，这样您就可以跟踪对于任何整数 n>=1，n*6+1 和 n*6+5 的素数。这可以节省三分之二的存储空间 space，但需要更多的簿记工作。您似乎试图从困难的版本开始。此外，如果涉及 Python 解释器评估 if 语句并进行列表的内存管理，所有簿记往往会变得昂贵。 Python/numpy 的数组切片是一种更自然的方法，而且对于您的目的来说它可能足够快。

关于您的问题：

(n%6==2) 是布尔值，将被解释为 0 或 1，并将再添加一个元素以防止出现 "off-by-one" 错误。
int(n**0.5)//3+1 应读作 int(int(n**0.5)/3) + 1。除法先于加法。除以三是为了让您只为 6n+1 和 6n+5 值分配 space。
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, ...].
将筛子元素设置为 False 是针对代表 6n+1 和 6n+5 的元素分别完成的，这就是为什么有两个赋值。
如果你在Python 2.7中是运行这个，整数除法总是被截断，所以9/2将得到4。在Python 3中，它会更好使用 // 运算符进行整数除法，即 (k*k)/3 应该是 (k*k)//3。

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

有人可以帮助我理解这个埃拉托色尼筛法脚本吗？最后几行让我难过

Could someone help me understand this Sieve of Eratosthenes script? The last couple lines have me stumped

python

primes

python-3.x