在Python中使用Atkin实现的Sieve

Question

我正在尝试在Wikipedia Link中实现Atkin的Sieve算法,如下所示:

阿特金筛子

到目前为止我尝试过的是Python中的代码:

import math
is_prime = list()
limit = 100
for i in range(5,limit):
    is_prime.append(False)

for x in range(1,int(math.sqrt(limit))+1):
    for y in range(1,int(math.sqrt(limit))+1):
        n = 4*x**2 + y**2

        if n<=limit and (n%12==1 or n%12==5):
            # print "1st if"
            is_prime[n] = not is_prime[n]
        n = 3*x**2+y**2
        if n<= limit and n%12==7:
            # print "Second if"
            is_prime[n] = not is_prime[n]
        n = 3*x**2 - y**2
        if x>y and n<=limit and n%12==11:
            # print "third if"
            is_prime[n] = not is_prime[n]

for n in range(5,int(math.sqrt(limit))):
    if is_prime[n]:
        for k in range(n**2,limit+1,n**2):
            is_prime[k] = False
print 2,3
for n in range(5,limit):
    if is_prime[n]: print n

现在我得到了错误

is_prime[n] = not is_prime[n]
IndexError: list index out of range

这意味着我正在访问列表中的值,其中索引大于List的长度.考虑x,y = 100时的条件,然后当然条件n = 4x ^ 2 + y ^ 2将给出大于列表长度的值.我在这里做错了吗？请帮忙.

编辑1 正如Gabe所建议的,使用

is_prime = [False] * (limit + 1)

注意:

for i in range(5,limit):
    is_prime.append(False)

确实解决了这个问题.

Answer 1

这是一个解决方案

import math

def sieveOfAtkin(limit):
    P = [2,3]
    sieve=[False]*(limit+1)
    for x in range(1,int(math.sqrt(limit))+1):
        for y in range(1,int(math.sqrt(limit))+1):
            n = 4*x**2 + y**2
            if n<=limit and (n%12==1 or n%12==5) : sieve[n] = not sieve[n]
            n = 3*x**2+y**2
            if n<= limit and n%12==7 : sieve[n] = not sieve[n]
            n = 3*x**2 - y**2
            if x>y and n<=limit and n%12==11 : sieve[n] = not sieve[n]
    for x in range(5,int(math.sqrt(limit))):
        if sieve[x]:
            for y in range(x**2,limit+1,x**2):
                sieve[y] = False
    for p in range(5,limit):
        if sieve[p] : P.append(p)
    return P

print sieveOfAtkin(100)

Answer 2

您的问题是您的限制是 100，但由于使用 初始化，您的is_prime列表中只有元素。limit-5range(5, limit)

由于此代码假设它可以访问最多limit索引，因此您需要limit+1在其中包含元素：is_prime = [False] * (limit + 1)

4x^2+y^2请注意，大于并不重要，limit因为它始终检查n <= limit。