# Python 3: Why is my prime sieve so slow?

I've implemented the sieve of Eratosthenes in Python 3 as follows:

def sieve_of_eratosthenes(n):
is_prime = [True] * (n+1)
is_prime[0] = False
is_prime[1] = False
p = 0

while True:
for i, prime in enumerate(is_prime):
if prime and i > p:
p = i
break
else:
break

multiple = p + p

while multiple <= n:
is_prime[multiple]= False
multiple = multiple + p

r = []
for i, prime in enumerate(is_prime):
if prime: r.append(i)
return r


Running this to 100,000 takes ~25 seconds.

I felt that was on the slow side, so I decided to take a different approach. The Sieve works prospectivly: it works ahead to clear out all multiples. I made a retrospective function, trial-dividing everything lower than the current prime-candidate:

def retrospect(n):
primes = [2, 3]
i = 5

while i <= n:
isprime = True
lim = math.sqrt(n)+1
for p in primes:
if p > lim:
break
if i % p == 0:
isprime = False
break
if isprime:
primes.append(i)
i += 2

return primes


This is way faster!

Questions:

• Are there any obvious shortcomings to the Sieve implementation?
• Is there a point where the retrospect gets slower than the sieve?

    while True:
for i, prime in enumerate(is_prime):
if prime and i > p:
p = i
break
else:
break


This is the cause of the slowness: track the state between loops and the 25 seconds drop to 0.05 seconds.

        multiple = p + p

while multiple <= n:
is_prime[multiple]= False
multiple = multiple + p


It's probably more Pythonic to use range here:

        for multiple in range(p + p, n + 1, p):
is_prime[multiple] = False


When handling larger inputs, it becomes worthwhile starting at p * p instead of p + p (on the basis that every composite multiple of p smaller than p * p also has a smaller prime factor).

    r = []
for i, prime in enumerate(is_prime):
if prime: r.append(i)
return r


IMO it makes more sense to inline this into the main loop.

Edited code at this point:

def sieve_of_eratosthenes(n):
is_prime = [True] * (n+1)
is_prime[0] = False
is_prime[1] = False

primes = []

for p in range(2, n+1):
if is_prime[p]:
primes.append(p)

for multiple in range(p * p, n + 1, p):
is_prime[multiple] = False

return primes


Now if you want more speed at the cost of complicating things ever so slightly, you can apply a simple wheel by special-casing the only even prime:

def sieve_of_eratosthenes(n):
is_prime = [True] * (n+1)
primes = [2]

for p in range(3, n + 1, 2):
if is_prime[p]:
primes.append(p)

for multiple in range(p * p, n + 1, p + p):
is_prime[multiple] = False

return primes

• Basically, your for loop only selects every number in range once and discards it if it isn't prime. The total number of iterations to find the next prime when n=1000 is 831, which when added to the number of primes < n is 999. In my original function, every number is tested again and again, leading to 77296 iterations. – steenbergh Jul 19 '18 at 11:04

Your prime sieve is overly complicated. You only need a single loop (and can also make it a generator to potentially save some memory):

def prime_sieve(limit):
prime = [True] * limit
prime[0] = prime[1] = False

for p, is_prime in enumerate(prime):
if is_prime:
yield p
for n in range(p * p, limit, p):
prime[n] = False


This also uses the fact that p + p is already known to be composite for most numbers (since it is 2 * p and all multiples of 2 have already been marked composite once you have added 2 to the primes) and indeed all multiples of n * p for n < p have already been excluded. So it is sufficient to start at p * p, because that is the first multiple that has not been marked of as composite.

This is strictly faster than your retrospect function:

It is also linear with n, whereas retrospect is roughly $\mathcal{O}(n^2)$, due to having to iterate over all already found primes.

In your other function I would also use a for loop. It has the nice feature that it has an else clause that gets executed if the loop was not interrupted with a break statement:

def retrospect(n):
primes = [2, 3]
lim = math.sqrt(n) + 1

for i in range(5, n + 1, 2):
for p in primes:
if p > lim:
break
if i % p == 0:
break
else:
primes.append(i)
return primes


It is also sufficient to calculate lim only once.

• Couldn't you take your last 2 lines and combine them into a slice assignment? prime[p*p:limit:p] = False should speed it up even more. – Sunny Patel Jul 19 '18 at 21:17
• @SunnyPatel That works only with advanced broadcasting rules like in numpy. Vanilla Python needs an iterable of correct length on the right side for that to work. – Graipher Jul 19 '18 at 21:20
• Oooh, I thought I had implemented one like this before in the same fashion. – Sunny Patel Jul 19 '18 at 21:22

Your sieve function is rather convoluted and slow. If we think carefully about what's actually required, for each identified prime number, we cross off multiples. We can stop when we get to $\sqrt{n}$ and each prime cross-out operation for a prime $i$ can start with $i^2$.

def sieve_of_eratosthenes(n):
is_prime = [True] * (n+1)
is_prime[0] = False
is_prime[1] = False

limit = int(math.sqrt(n))+1
for i in range(2, limit):
if is_prime[i]:
for j in range(i*i, n+1, i):
is_prime[j] = False

r = []
for i, prime in enumerate(is_prime):
if prime: r.append(i)
return r


A still faster way to do it would be to note that all even numbers except for 2 are composite (that is, non-prime). We can slightly tweak the code to account for that:

def sieve_of_eratosthenes(n):
is_prime = [False, True] * int(n/2 + 1)
del is_prime[n+1:]
is_prime[1] = False
is_prime[2] = True

limit = int(math.sqrt(n))+1
for i in range(3, limit):
if is_prime[i]:
for j in range(i*i, n+1, 2*i):
is_prime[j] = False

r = []
for i, prime in enumerate(is_prime):
if prime:
r.append(i)

return r


## Performance

Just for fun, I decided to repeat and expand the performance measurement that @Graipher made. Here's the code I used to do it:

import time
import matplotlib.pyplot as plt

def time_algo(algo, times):
result = []
for t in times:
start = time.time()
algo(t)
result.append((time.time() - start)*1000)

return result

names = [
#'retrospect',
taylor1,
graipher,
edward1,
taylor2,
edward2,
]
numarray = [1000, 3000,
10000, 30000,
100000, 300000,
1000000, 3000000,
10000000, 30000000,]
trials = []
for algo in names:
trials.append(time_algo(algo, numarray))
for t in trials:
plt.plot(numarray, t, '.-')
plt.ylabel('time (ms)')
plt.xlabel('n')
plt.title("Prime algorithm perfomance")
plt.legend([alg.__name__ for alg in names])
plt.show()


Here are the results from the test on my machine:

• I've corrected the limit, as noted by @ack. – Edward Jul 19 '18 at 14:54

@Edward: corrected limit

def sieve_of_eratosthenes(n):
is_prime = [False, True] * int(n/2 + 1)
n_plus_one = n + 1  # save n+1 for later use
del is_prime[n_plus_one:]  #
is_prime[1] = False
is_prime[2] = True

#limit = int(math.sqrt(n))
limit = int(n**.5) + 1  # import math  not required, limit corrected!
for i in range(3, limit):
if is_prime[i]:
for j in range(i*i, n_plus_one, 2*i):  # limit corrected!
is_prime[j] = False

r = []
for i, prime in enumerate(is_prime):
if prime:
r.append(i)

return r


without corrected limit s_o_e(9) yields [2,3,5,7,9], so 9 would be prime (as would 11,13,17,...)

• I don't think this constitutes an answer, it should probably be a comment on @Edward answer. Once you have enough reputation (50), you will be able to comment everywhere. – Graipher Jul 19 '18 at 7:36
• @Graipher I agree, though this is a valid point; in fact I believe ack deserves rep for this find. It's a shame to see the user not rewarded for their effort because of a lack of reputation-points... – steenbergh Jul 19 '18 at 7:56
• @steenbergh I agree that it is a good find. But if they had had more rep and commented this instead, they would not have gotten any rep from upvotes on that comment either. – Graipher Jul 19 '18 at 7:58
• @Graipher I think this is a significant improvement enough to Edward's answer that is possible to stand as its own answer. (For this time) – Simon Forsberg Jul 19 '18 at 8:07
• I suspect there is little difference among those options, but the way to find out is to measure. – Edward Jul 20 '18 at 23:08

On the performance side, using numpy arrays makes the computations 10 times faster of course.

Using Edward's code for performance, together with the following code ("Edward1" algorithm adapted to numpy) :

def sieve_of_eratosthenes(n):
import numpy as np

is_prime = np.ones(n,dtype=bool)
is_prime[0:2] = False

limit = int(math.sqrt(n))+1
for i in range(2, limit):
if is_prime[i]:
is_prime[np.arange(i*i,n,i)] = False

r = np.where(is_prime == True)[0]
return r


I get on my machine :

• Will definitely look into Numpy then! – steenbergh Jul 19 '18 at 13:39
• This reproduces the bug in Edward's code. It's also not actually reviewing the code in the question, so by the standards of this site it shouldn't be an answer. – Peter Taylor Jul 19 '18 at 13:40
• Sorry about that. This has just been corrected for. – D. Le Borgne Jul 19 '18 at 13:41