Understanding Computational Efficiency By Writing a Prime Number Sieve

Feb 21, 2021 · 9 min read · Python Prime Numbers Math

UMAP projection of the prime numbers. — What do prime numbers 'look' like?

One of the earlier exercises in computational methods is the prime number sieve. The straightforward yet non-intuitive nature of prime numbers makes it a mathematical treat. At first, the concept is simple:

$$ \begin{array}{l} \text{for} \quad p, q \in \mathbb{N} \\ \text{s.t.} \quad p \geq 2 \quad \text{and} \quad q \neq p,1 \\ p \quad \text{is prime} \iff \frac{p}{q} \notin \mathbb{N} \end{array} $$

This may not be the most graceful notation, but I will lay it out more clearly. Given $p$ and $q$ in $\mathbb N$ where $\mathbb N$ is the set of 'natural numbers $\mathbb N = \{0,1,2,3,...\} $. Such that $p$ is greater than or equal to 2, and $q$ is not $p$ or $1$. $p$ is prime if and only if $p$ divided by $q$ is not a natural number. The only divisors of a prime number are $p$ and $1$.

$$ \lim_{n\rightarrow \infty} \frac{\pi(n)}{\frac{n}{\ln{n}}} = 1 $$

Where $\pi(n)$ is the number of primes $p$ such that $2\leq p \leq n$. This can be used to assist in testing if a number is prime, but out of scope for the content of this article.

Regardless of the tricks that have been developed, I still won't be able to help a four year old solve if $10,001$ is prime or not. That is what motivates me to create these simple computer programs. First you can solve answer the question, but equally investigating on how to improve the efficiency of a program also helps refine your understanding of the problems and ultimately create better code/math as a result.

Let's write a prime number sieve.

Pseudocode:

 1set upper limit
 2make a prime list
 3
 4loop over n from 2 to upper limit
 5    
 6    set "still_prime" flag to True
 7    
 8    loop on p from 2 to n-1
 9            
10        is n divisible by p?
11            yes: set "still_prime" flag to False
12    is "still_prime" flag True?
13        yes: append n to prime list
14output primes found

This is the most intuitive understanding of the concept, how about implementing it in Python. Note I'm not going to use Numpy, although I could squeeze some improvements in, but maybe lose some in other places as well! This program only has time for measurement and uses % the modulus function for efficient computing.

Simple Prime Number Sieve

 1import time
 2
 3upper_limit = 10001
 4primes_found = []
 5
 6time_start = time.time()
 7# iterate from 2 to limit
 8for n in range(2,upper_limit+1):
 9
10    still_prime = True
11    # iterate from 2 to number
12    for p in range(2,n):
13        # if the number is divisible, not prime.
14        if n%p == 0:
15            
16            still_prime = False
17
18    if still_prime:
19        
20        primes_found.append(n)
21        
22time_end = time.time()
23
24print( "Total number of primes found:", len(primes_found) )
25print( "Primes found:" )
26#print last 10 primes
27print( primes_found[-10:] )
28print("Run time:", time_end-time_start)

1Total number of primes found: 1229
2Primes found:
3[9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
4Run time: 5.7099528312683105

It works, it is correct, but at the cost of 5.71 seconds on my machine! Why invest in a nice gaming computer in the first place?! Sheesh. I got problems to solve on Project Euler and I am sure as hell not going to wait 5.71 seconds just to find that read the problem wrong. We. can. do. better.

We have two for loops iterating across large numbers... the statement if n%p = 0: still_prime = False really only needs to be called if the 'still_prime' condition is still true.

Modification #1: Check if previous number was prime before testing

 1upper_limit = 10001
 2primes_found = []
 3
 4time_start = time.time()
 5
 6for n in range(2,upper_limit+1):
 7
 8    still_prime = True
 9
10    for p in range(2,n):
11        #Added this condition
12        if still_prime:
13            
14            if n%p == 0:
15                
16                still_prime = False
17
18    if still_prime:
19        
20        primes_found.append(n)
21
22
23time_end = time.time()
24
25print( "Total number of primes found:", len(primes_found) )
26print( "Primes found:" )
27#print last 10 primes
28print( primes_found[-10:] )
29print("Run time:", time_end-time_start)

1    Total number of primes found: 1229
2    Primes found:
3    [9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
4    Run time: 3.2463064193725586

A substaintial reduction of computing time, from 5.70 to 3.25 seconds, while continuing to perform correctly. But upon further investigation, there is more room for improvement. In fact, the most apparent truth of primes is regarding what they can be composed of. Let's say $n$ is 100. 100 = $10*10$. That means that if I want to see if a number is prime less than 100, I only have to check to the square root of 100, 10!

Brilliant.

Modification #2: Teset only to $\sqrt n$

 1upper_limit = 10001
 2primes_found = []
 3
 4time_start = time.time()
 5
 6for n in range(2,upper_limit+1):
 7
 8    # set sqrt(n)
 9    sqrt_of_n = n**0.5
10    # set "checked all" flag to False
11    checked_all = False
12
13    still_prime = True
14
15    for p in range(2,n):
16        # add the 'not checked_all' to the condition statement
17        if still_prime and not checked_all:
18            
19            if p > sqrt_of_n:
20                # check if you have covered the numbers to sqrt(n) 
21                checked_all = True
22                
23            elif n%p == 0:
24                
25                still_prime = False
26
27    if still_prime:
28        primes_found.append(n)
29
30
31time_end = time.time()
32
33print( "Total number of primes found:", len(primes_found) )
34print( "Primes found:" )
35#print last 10 primes
36print( primes_found[-10:] )
37print("Run time:", time_end-time_start)

1Total number of primes found: 1229
2Primes found:
3[9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
4Run time: 2.696723461151123

Next, now that we are computing to $\sqrt{n}$, we can also realize that we don't need to check p from 2 to n. We only need to check to the primes found for our second for loop.

Modification #3: Test from current primes_list, not to n.

 1upper_limit = 10001
 2primes_found = []
 3
 4time_start = time.time()
 5
 6for n in range(2,upper_limit+1):
 7
 8    sqrt_of_n = n**0.5
 9    checked_all = False
10    still_prime = True
11
12    #changed range from 2,n to primes_found
13    for p in primes_found:
14
15        if still_prime and not checked_all:
16            
17            if p > sqrt_of_n:
18
19                checked_all = True
20            
21            elif n%p == 0:
22            
23                still_prime = False
24
25    
26    if still_prime:
27    
28        primes_found.append(n)
29
30time_end = time.time()
31print( "Total number of primes found:", len(primes_found) )
32print( "Primes found:" )
33#print last 10 primes
34print( primes_found[-10:] )
35print("Run time:", time_end-time_start)

1Total number of primes found: 1229
2Primes found:
3[9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
4Run time: 0.25299906730651855

The time reduced from 2.80 to 0.25 seconds just by iterating over the list and not n. That is over an order of magnitude improvement! Now the last major improvement... While the inside for loop change has dramatically sped up the code, we are still continuing to iterate over numbers when the 'still_prime' and 'checked_all' flags are checked... This isn't efficient as we have already met our condition statements once those happen. We should exit the loop after we have searched all possible choices.

Modification #4 Exit loop if non-prime is found or no more divisors exist.

 1upper_limit = 10001
 2primes_found = []
 3
 4time_start = time.time()
 5
 6for n in range(2,upper_limit+1):
 7
 8    sqrt_of_n = n**0.5
 9    checked_all = False
10    still_prime = True
11
12    for p in primes_found:
13
14        if still_prime and not checked_all:
15
16            if p > sqrt_of_n:
17            
18                checked_all = True
19            elif n%p == 0:
20                still_prime = False
21        
22        # else!
23        else:
24            # exit from the loop
25            break
26
27    if still_prime:
28    
29        primes_found.append(n)
30
31time_end = time.time()
32
33print( "Total number of primes found:", len(primes_found) )
34print( "Primes found:" )
35#print last 10 primes
36print( primes_found[-10:] )
37print("Run time:", time_end-time_start)

1Total number of primes found: 1229
2Primes found:
3[9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
4Run time: 0.01798725128173828

Amazing, another order of magnitude of computation time reduction all because of exiting a loop. Recall that the original algorithm time was 5.71 seconds, now reduced to 0.018 seconds. This is done on my semi-poweful desktop, and other individuals might have different results, but the improvements are real!

This excercise makes it very clear to see how even the simplest of computer code can quickly turn ugly and inefficient! All through the better understanding of the mathematics behind the problem and the engineering techniques to help execute.

**the internet has entered the chat* *

While the method that we have developed is certainly fast and we have made some great improvements, it is not the fastest by any means. I think the more valuable lesson is the subtle modification of code to extract over two orders of magnitude of performance increase! Upon finalizing this post, I decided to see what the rest of the internet community had to share on 'fastest prime sieves'. I was not disappointed when I came across this stackoverflow post. Now the input limits are constrained to N>=6 the succinctness and application of very non-intuitive instructions makes this code lightning fast. This is the code I am using now for all of my Project Euler problems now!

 1def primesbelow(N):
 2    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
 3    #""" Input N>=6, Returns a list of primes, 2 <= p < N """
 4    correction = N % 6 > 1
 5    N = {0:N, 1:N-1, 2:N+4, 3:N+3, 4:N+2, 5:N+1}[N%6]
 6    sieve = [True] * (N // 3)
 7    sieve[0] = False
 8    for i in range(int(N ** .5) // 3 + 1):
 9        if sieve[i]:
10            k = (3 * i + 1) | 1
11            sieve[k*k // 3::2*k] = [False] * ((N//6 - (k*k)//6 - 1)//k + 1)
12            sieve[(k*k + 4*k - 2*k*(i%2)) // 3::2*k] = [False] * ((N // 6 - (k*k + 4*k - 2*k*(i%2))//6 - 1) // k + 1)
13    return [2, 3] + [(3 * i + 1) | 1 for i in range(1, N//3 - correction) if sieve[i]]

 1upper_limit = 10001
 2primes_found = []
 3
 4
 5time_start = time.time()
 6primes_found = primesbelow(upper_limit)
 7time_end = time.time()
 8
 9print( "Total number of primes found:", len(primes_found) )
10print( "Primes found:" )
11#print last 10 primes
12print( primes_found[-10:] )
13print("Run time:", time_end-time_start)

1Total number of primes found: 1229
2Primes found:
3[9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973]
4Run time: 0.0010001659393310547

0.001 seconds for 10001 numbers. Another order of magnitude reduced due to some very clever arrangement of what is important in the prime number problem. This is the great value of the community of nerds. It is worth mentioning that none of the algorithms are parallel processing, quantum computing, next-gen architecture. But we can still continue to search for the most efficient methods and hopefully apply them to daily life in meaningful ways. For example, the UMAP method that created the feature image of this post where the prime number list is projected to a 2D space to show some fascinating geometrical structures that have yet to be fully understood!

As always, I'm hoping to improve on this work. Thanks for reading!

Consider this problem solved.

QED