12
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – John Davis Jun 11 '14 at 15:38
  • \$\begingroup\$ 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 \$\endgroup\$ – sundar nataraj Jun 12 '14 at 5:56
  • \$\begingroup\$ i want to compare speed of both ur code and this one codereview.stackexchange.com/questions/6477/… \$\endgroup\$ – sundar nataraj Jun 12 '14 at 6:32
3
\$\begingroup\$
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.

\$\endgroup\$
  • \$\begingroup\$ 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! \$\endgroup\$ – BeetDemGuise Jun 12 '14 at 12:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.