# Optimise Sieve of Eratosthenes

I have written the following code for finding out the prime numbers less than or equal to n. When I enter the input n = 10^7, my program crashes. So I was wondering whether my code can be further optimised.

def sieve(n):
array = []
for i in range(1, n):
array.append(True)
for i in range(0, n-1):
for j in range(2*i+4, n+1, i+2):
if array[j-2] != False:
if j % (i+2) == 0:
array[j-2] = False
final = []
for i in range(0, n-1):
if array[i] == True:
final.append(i+2)
return final


First, I would clean up the indices in the inner loop. I think the +2 and -2 are a bit confusing at the first glance.

def sieve(n):
array = []
for i in range(0, n + 1):
array.append(True)
for i in range(2, n + 1):
for j in range(2*i, n+1, i):
if array[j] != False:
if j % (i) == 0:
array[j] = False
final = []
for i in range(0, n+1):
if array[i] == True:
final.append(i)
return final


This way, array indices directly correspond to their number, no conversion needed. This wastes array space for 0 and 1, but I think it is a lot easier to understand.

The primary optimization that is possible adjusts the ranges of the loops. Each non-prime <= n has a divisor <= sqrt(n) and thus we can limit the outer loop to range(2, math.sqrt(100) + 1)

Additionally we don't actually need to run the inner loop for numbers which are non-prime (since each non-prime has prime divisors) and we can add an if array[i] == True: to reduce the number of inner loops further.

The range of the inner loop can also be reduced. It actually can start at i^2 instead of 2*i since the argument made earlier also applies here. All non-primes smaller than i^2 must have an divisor < i and thus were already set to false in an earlier iteration of the outer loop.

If we apply these changes, we get the following code:

def faster_sieve(n):
array = []
for i in range(0, n + 1):
array.append(True)

for i in range(2, int(math.sqrt(n)) + 1):
if array[i] == True:
for j in range(i*i, n + 1, i):
if array[j] != False:
if j % i == 0:
array[j] = False
final = []
for i in range(2, n + 1):
if array[i] == True:
final.append(i)
return final


For comparison we can run

faster_sieve(int(math.pow(10, 7)))


and

sieve(int(math.pow(10, 7)))


On my machine:

~ time python faster_sieve.py
python faster_sieve.py  5.44s user 0.04s system 99% cpu 5.477 total

~ time python sieve.py
python sieve.py  32.97s user 0.03s system 99% cpu 33.003 total


Which is a lot faster!

Now we could to some python microoptimizations, but I don't know anything about python. A first step could be to replace the append-loop with something faster like array = [True] * (n + 1), which saves another second on my machine

~ time python faster_sieve.py
python faster_sieve.py  4.46s user 0.02s system 99% cpu 4.475 total


So, yes indeed, your code could be further optimized.

/e I might add that there are already good reviews of python code for the sieve on this site. For example this which applies a lot of python optimization.

/e2

Looking at the code again and looking at Wikipedia, the checks inside the inner loop are actually not needed since j is per loop definition a multiple of i and we can simplify it to

for j in range(i*i, n + 1, i):
array[j] = False


which optimizes the program further:

~ time python faster_sieve.py
python faster_sieve.py  2.61s user 0.01s system 99% cpu 2.617 total


/e3 After reading jerry's response and finding two bugs in his implementation, I think there is another easy optimization that can be added to his approach. Since i is always odd in his inner loop, so is i * i. Thus we can increase the increment step of the inner loop from i to i * 2 (i * i + i would be even). This result in the following results:

Simon: 3.5496606826782227
Jerry (broken): 1.7031567096710205
Josay: 2.2623558044433594
Krypton: 5.4344189167022705
Krypton2: 2.5819575786590576
Jerry2: 1.396036148071289


Krypton2 and Jerry2 being:

def krypton2(n):
array = [True] * (n + 1)

for i in range(2, int(math.sqrt(n)) + 1):
if array[i] == True:
for j in range(i*i, n + 1, i):
array[j] = False
final = []
for i in range(2, n + 1):
if array[i] == True:
final.append(i)
return final

def jerry2(n):
array = [True] * n

result = []

result.append(2)
for i in range(4, n, 2):
array[i] = False;

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

for i in range(3, n, 2):
if array[i]:
result.append(i)
return result

• Nice additions. I hadn't really considered i*i +i being even--excellent observation. I estimate that one more week at the rate we're going, and we'll reach the point of: "a single properly timed NOP produces all the primes from 1 to N". ;-) – Jerry Coffin Oct 6 '16 at 0:10

Initialisation of the list

array = []
for i in range(1, n):
array.append(True)


calls append many times. There might be a better way. You could use list comprehension like this:

array = [True for _ in range(1, n)]


but you could also use the * operator:

array = [True] * (n-1)


Boolean checks

You could rewrite: if array[j-2] != False as if array[j-2] and if array[i] == True as if array.

Variable name

array is a pretty bad name for an array. It could be a good idea to give a name that tells what's interesting about this array's content. As you use it to know whether a value is prime, I'd call it is_prime but you might find a better name.

List comprehension (again)

You could rewrite the final section of your function using a list comprehension.

    return [i + 2 for i in range(0, n-1) if is_prime[i]]


Loop like a native

When you use the index to iterate over an array in Python, there is usually a better way to do what you want to achieve. I highly recommend Ned Batchelder's presentation called : Loop Like A Native to know more about this.

In your case, you could write:

return [i + 2 for i, p in enumerate(is_prime) if p]


Better algorithm

The whole point of the sieve is not to use division (or the modulo operator).

I'd use something like:

import math

def sieve(n):
is_prime = [True] * n
is_prime[0] = is_prime[1] = False
for i in range(2, int(math.sqrt(n)) + 1):
if is_prime[i]:
for j in range(i * i, n, i):
is_prime[j] = False
return [i for i, p in enumerate(is_prime) if p]

• > calls append many times. There might be a better way. You could use list comprehension like this Are you sure list comprehension does not use append under the hood though? – Daerdemandt Oct 5 '16 at 16:53
• I'm pretty sure it's faster (it calls something like append but does not need to perform the lookup multiple times) but haven't performed a benchmark to prove so. I'm currently reading stackoverflow.com/questions/22108488/… which should answer the question. – SylvainD Oct 5 '16 at 16:55
• Yeah, I remember reading that exact thread. I'd expect list comprehension to call append exactly the same number of times as a for loop, but be marginally faster because of reasons listed in thread you linked - whereas array multiplication would indeed be way faster than either of those. – Daerdemandt Oct 5 '16 at 17:10

You've gotten a number of suggestions about how to optimize the code already, but I can't say any of them leaves me feeling like the code is really as clean or simple as it could be. I haven't spent a lot of time on trying to profile to quantify the precise effect on execution speed, but it seems to me that it's pretty straightforward to do the job something like this:

def sieve(n):
array = [True] * n

result = []

result.append(2)
for i in range(4, n, 2):
array[i] = False;

for i in range(3, int(math.sqrt(n)), 2):
if (array[i]):
result.append(i)
for j in range (i*i, n, i):
array[j] = False;

for i in range(int(math.sqrt(n)), n):
if (array[i]):
result.append(i)


This special-cases 2, because it's, well, a special case: the only even prime number (so special, that lots of papers about primes and their properties only mention 2 by exclusion (e.g., "let N be an odd prime"). So, we start by adding 2 to the list of primes, and crossing out all the multiples of 2 (if we wanted to get tricky, we wouldn't even allocate space in array for multiples of 2 at all).

After that, we don't need to consider even numbers at all, so we start from 3, and only generate odd numbers as candidates.

We only need to consider numbers up to the square root of our limit, because any factor larger than the square root has to match up with another factor smaller than the square root.

For each prime we find, we only need to cross out numbers starting from it square because multiples of all numbers smaller than it have already been crossed out (e.g., 2N would already have been crossed out as a multiple of 2, 3N as a multiple of 3, and so on).

To avoid re-traversing the result list after we've found the primes, we just add each to the output as we find them. Not sure if this is really any faster, but it avoids arguments about how to do list traversal.

A quick benchmark indicates that this at least slightly faster than any of the others I see posted at the present time:

import time
import math

def jerry(n):
array = [True] * n

result = []

result.append(2)
for i in range(4, n, 2):
array[i] = False;

for i in range(3, int(math.sqrt(n)), 2):
if (array[i]):
result.append(i)
for j in range (i*i, n, i):
array[j] = False;

for i in range(int(math.sqrt(n)), n):
if (array[i]):
result.append(i)

def sieve(li, x):
for i in range(2*x, len(li), x):
li[i] = False
return li

def simon(limit):
li = [False, True]*(limit//2)
li[1] = False
li[2] = True
for i in range(3, limit, 2):
if li[i]:
li = sieve(li, i)
primes = [x for x in range(len(li)) if li[x]]
return primes

def josay(n):
is_prime = [True] * n
is_prime[0] = is_prime[1] = False
for i in range(2, int(math.sqrt(n)) + 1):
if is_prime[i]:
for j in range(i * i, n, i):
is_prime[j] = False
return [i for i, p in enumerate(is_prime) if p]

def krypton(n):
array = []
for i in range(0, n + 1):
array.append(True)

for i in range(2, int(math.sqrt(n)) + 1):
if array[i] == True:
for j in range(i*i, n + 1, i):
array[j] = False
final = []
for i in range(2, n + 1):
if array[i] == True:
final.append(i)
return final

if __name__ == '__main__':
n = 10000000

tests = [simon, jerry, josay, krypton]
names = ["Simon", "Jerry", "Josay", "Krypton"]

for i in range(0, len(tests)):
start = time.time()
result = tests[i](n)
end = time.time()
print '{}: {}'. format(names[i], (end-start))


Results:

Simon: 3.54199981689
Jerry: 2.21599984169
Josay: 3.52499985695
Krypton: 5.94400000572


Note: after posting these, @Krypt0n has added a couple more possibilities that look quite good but aren't (at least yet) included here.

• Could you update my code with the two changes that are described in my answer after the full code snippet? The performance should not be that much worse than the other solutions. – Krypt0n Oct 5 '16 at 22:16
• Does your solution add anything above sqrt(n) to the result list? – Krypt0n Oct 5 '16 at 22:27
• @Krypt0n: I'll work on adding those changes when I get a chance. And....oops. I need to work on my code a bit too... – Jerry Coffin Oct 5 '16 at 22:41

I have some background in math so a further improvement would be

def sieve(li, x):
for i in range(2*x, len(li), x):
li[i] = False
return li

def get_primes(limit):
li = [False, True]*(limit//2)
li[1] = False
li[2] = True
for i in range(3, limit, 2):
if li[i]:
li = sieve(li, i)
primes = [x for x in range(len(li)) if li[x]]
return primes


because no even number beside 2 is prime. It probable is faster in one function. And will be even faster as a function in i class. Don't know why class methods is faster then functions..

Can't explain this either, but it's faster to iterate over range objects then lists of numbers. I'm referring to the second last row in the code. But iterating over a list of numbers is faster then enumeration.

When doing the first 50 project euler problem, I was keen on doing them in less then 1 second. This kind of things add up.