# Sieve of Atkins algorithm

I got the inspiration to write this algorithm from this question. The OP cited the Sieve of Atkin as one of the fastest ways to generate primes. So I figured I would give it a try:

def sieve_primes(limit=100):
sieve = [False]*limit
loop_range = range(1, int(limit**(1/2)))

for x, y in [(x,y) for x in loop_range for y in loop_range]:
for x_coef, y_coef, mods in [(4, 1, [1, 5]), (3, 1, ), (3, -1, )]:
if y_coef == -1 and x <= y:
continue

value = x_coef * x**2 + y_coef * y**2
if value > limit:
continue

if value % 12 in mods:
sieve[value] = not sieve[value]

for value in loop_range[4:]:
if sieve[value]:
for square in [i * value**2 for i in range(limit//value**2+1)
if i * value**2 < limit]:
sieve[square] = False

return sieve


The algorithm is based off of pseudocode given by the reference to the Sieve of Atkins above. The reference does say:

This pseudocode is written for clarity. Repeated and wasteful calculations mean that it would run slower than the sieve of Eratosthenes. To improve its efficiency, faster methods must be used to find solutions to the three quadratics.

I wanted to make the pseudocode feel as Pythonic as possible. I also tried to optimize it as much as I could without delving too much into the mathematics.

Do you guys see any improvements, performance or style-wise?

Here's some profiling:

import timeit

for x in [100, 1000, 10000]:
time = timeit.timeit('sieve_primes({})'.format(x),
setup='from __main__ import sieve_primes', number=10000)
print('Number of primes found: {}'.format(x))
print('\tTotal time: {}\n\tTime per iteration: {}'.format(time, float(time)/10000))

# Results
Number of primes found: 100
Total Time: 2.3724100288364544
Time per Iteration: 0.00023724100288364545
Number of primes found: 1000
Total Time: 26.933066114727716
Time per Iteration: 0.0026933066114727716
Number of primes found: 10000
Total Time: 297.5023943142171
Time per Iteration: 0.0297502394314171

• You might find some interesting information on my question.. THere is a sieve of erasthones and atkins in the answers, and they got some wicked fast times. Its in C#, but you might still get something from it. Jun 11, 2014 at 15:38
• Wow! u have coded the algorithm ! Great Man thank you but still it is not producing primes i guess.how to produce primes using this code Jun 12, 2014 at 5:56
• i want to compare speed of both ur code and this one codereview.stackexchange.com/questions/6477/… Jun 12, 2014 at 6:32

for x, y in [(x,y) for x in loop_range for y in loop_range]:


Here you could avoid the list creation by using a generator (by changing the [] to ()) or you could use itertools.product. Probably no difference performance-wise, but the latter is less code.

    for x_coef, y_coef, mods in [(4, 1, [1, 5]), (3, 1, ), (3, -1, )]:


Again, you are creating lists which takes time. This whole thing you iterate over could be a constant nested tuple that's only allocated once at the top (or even outside the function).

        for square in [i * value**2 for i in range(limit//value**2+1)
if i * value**2 < limit]:


At this point you can probably guess that I'll recommend doing away with the list. A generator should work. If on Python 2, xrange is probably better than range as well.

(Note, I didn't mention loop_range. That could be changed into xrange invocations if you are on Python 2, but you seem to be on Python 3.)

There's also a question mark I have: int(limit**(1/2)) – that works in Python 3, but not Python 2, because 1/2 == 0. Since your code seemed to work for you, I assume you are on Python 3. I'm not sure if Python 3 people think making it portable with a dot/cast is a good idea or not (i.e. 1./2). There's always the explicit sqrt if you can't decide.

• You were right to assume I'm using Python 3. If this wasn't a 'for-fun' algorithm, I would be more concerned about backwards-compatibility. Thanks for the input! Jun 12, 2014 at 12:04