Sieve Of Atkin Implementation in Python

Question

I am trying to implement the algorithm of Sieve of Atkin given in Wikipedia Link as below:

Sieve Of Atkin

What I've tried so far is the implementation in Python given by following Code:

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

Now I get error as

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

this means that I am accessing the value in list where the index is greater than length of List. Consider the Condition when x,y = 100, then of-course the condition n=4x^2+y^2 will give value which is greater than length of list. Am I doing something wrong here? Please help.

EDIT 1 As suggested by Gabe, using

is_prime = [False] * (limit + 1)

insted of :

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

did solved the problem.

Answer 1

You problem is that your limit is 100, but your is_prime list only has limit-5 elements in it due to being initialized with range(5, limit).

Since this code assumes it can access up to limit index, you need to have limit+1 elements in it: is_prime = [False] * (limit + 1)

Note that it doesn't matter that 4x^2+y^2 is greater than limit because it always checks n <= limit.

Answer 2

Here is a solution

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 3

This is the optimized implementation proposed by Zsolt KOVACS:

    import math
    import sys
    
    def sieveOfAtkin(limit):
        P = [2,3]
        r = range(1,int(math.sqrt(limit))+1)
        sieve=[False]*(limit+1)
        for x in r:
            for y in r:
                xx=x*x
                yy=y*y
                xx3 = 3*xx
                n = 4*xx + yy
                if n<=limit and (n%12==1 or n%12==5) : sieve[n] = not sieve[n]
                n = xx3 + yy
                if n<=limit and n%12==7 : sieve[n] = not sieve[n]
                n = xx3 - yy
                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]:
                xx=x*x
                for y in range(xx,limit+1,xx):
                    sieve[y] = False
        for p in range(5,limit):
            if sieve[p] : P.append(p)
        return P
    
    primes = sieveOfAtkin(int(sys.argv[1]))    
    print (primes)

You pass upper limit as the first argument. This program runs in about 6s on my machine comparing to the original which runs in 21s for 10 million limit. What I did:

• replaced exponentiation with multiplication
• precalculated some multiplications