Why is my implementation of the Sieve of Atkin overlooking numbers close to the specified limit?

Question

My implementation of Sieve of Atkin either overlooks primes near the limit or composites near the limit. while some limits work and others don't. I'm am completely confused as to what is wrong.

def AtkinSieve (limit):
results = [2,3,5]
sieve = [False]*limit
factor = int(math.sqrt(lim))
for i in range(1,factor):
    for j in range(1, factor):
        n = 4*i**2+j**2
        if (n <= lim) and (n % 12 == 1 or n % 12 == 5):
            sieve[n] = not sieve[n]
        n = 3*i**2+j**2
        if (n <= lim) and (n % 12 == 7):
            sieve[n] = not sieve[n]
        if i>j:
            n = 3*i**2-j**2
            if (n <= lim) and (n % 12 == 11):
                sieve[n] = not sieve[n]
for index in range(5,factor):
    if sieve[index]:
        for jndex in range(index**2, limit, index**2):
            sieve[jndex] = False
for index in range(7,limit):
    if sieve[index]:
        results.append(index)
return results

For example, when I generate a primes to the limit of 1000, the Atkin sieve misses the prime 997, but includes the composite 965. But if I generate up the limit of 5000, the list it returns is completely correct.

Answer 1

• Change lim to limit. Of course you must have known that.
• Since sieve = [False]*limit, the largest index allowed is limit-1.

  However, on this line

  if (n <= limit) and (n % 12 == 1 or n % 12 == 5):
  

  you are checking if n<=limit. If n==limit then sieve[n] raises an IndexError. Try your algorithm with a small value of limit (e.g. n=50). You'll see this error come up. An easy fix is to use

  sieve = [False]*(limit+1)
  

  The easy fix is a bit wasteful since sieve[0] is never used. So you might think a better fix is to keep sieve = [False]*limit, but fix all your other code by stepping the index on sieve down by one. (E.g., change sieve[n] to sieve[n-1] everywhere, etc.) However, this will force you to do a number of extra subtractions which will not be good for speed. So the easy/wasteful solution is actually probably the better option.

• According to http://en.wikipedia.org/wiki/Sieve_of_Atkin, x should be an integer in [1,sqrt(limit)], inclusive of the endpoints.

  In your code

  factor = int(math.sqrt(limit))
  

  and int takes the floor of math.sqrt(limit). Furthermore,

  range(1,factor) goes from 1 to factor-1. So you are off by 1.

  So you need to change this to

  factor = int(math.sqrt(limit))+1
  

• See Fastest way to list all primes below N for an alternative (and faster) implementation of the Sieve of Atkin, due to Steve Krenzel.

def AtkinSieve (limit):
    results = [2,3,5]
    sieve = [False]*(limit+1)
    factor = int(math.sqrt(limit))+1
    for i in range(1,factor):
        for j in range(1, factor):
            n = 4*i**2+j**2
            if (n <= limit) and (n % 12 == 1 or n % 12 == 5):
                sieve[n] = not sieve[n]
            n = 3*i**2+j**2
            if (n <= limit) and (n % 12 == 7):
                sieve[n] = not sieve[n]
            if i>j:
                n = 3*i**2-j**2
                if (n <= limit) and (n % 12 == 11):
                    sieve[n] = not sieve[n]
    for index in range(5,factor):
        if sieve[index]:
            for jndex in range(index**2, limit, index**2):
                sieve[jndex] = False
    for index in range(7,limit):
        if sieve[index]:
            results.append(index)
    return results