trjhtr

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

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 $mathcalO(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.

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 $sqrtn$ 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:

@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,...)

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 :