Sieve of Atkins algorithm

Question

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, [7]), (3, -1, [11])]:
            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

Answer 1

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, [7]), (3, -1, [11])]:

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.