The Sieve of Eratosthenes
Nicky Meuleman
Posted on August 13, 2021
The sieve of Eratosthenes finds all prime numbers up to a given limit.
Method
The algorithm starts out by assuming all numbers are prime, and marking them as such.
At the end of the algorithm, only prime numbers up to an upper limit will still be marked.
The number 1 is a special case, so we start off by unmarking it.
Then we go through the numbers one by one.
For every non-prime number we find, skip to the next number.
If a number is still marked as prime when we get to it, that means it is prime.
Before moving on to the next number, we first unmark every multiple of the found prime.
Those multiples can be divided through the prime number we just found, so by definition isn't prime.
We repeat this process until we reach the upper limit.
Every number that is still marked as prime, is truly prime.
Optimizations
By using some math we can do significantly less work while still getting the same result.
Repeat until the square root
TL;DR: Only check numbers up to the square root of the upper limit.
After that, every number up to that limit will be accurately marked, because math is cool.
While iterating through all numbers, we can stop at the square root of the upper limit.
Any non-prime can be expressed as the product of 2 numbers that are not 1 or itself.
n = a * b
a and b are factors of n.
n = sqrt(n) * sqrt(n), so one factor has to be less than or equal to sqrt(n) while the other is greater than or equal to that square root.
a <= sqrt(n) <= b
Up to any number n, all multiples of a number bigger than sqrt(n) must have a factor smaller than sqrt(n).
As a result that multiple will already be unmarked.
The factors can both be identical to sqrt(n).
49 for example
This means that all the non-primes >= sqrt(limit) will be unmarked in the process of checking every number <= sqrt(limit).
Example
sqrt(21) = 4.58
Any number up to 21 that is a multiple of a number larger than 4.58 will have a factor smaller than 4.58.
Because 18 is a number up to 21.
It is also a multiple of a number that is bigger than 4.58.
That means a factor of 18 must be smaller than 4.58.
That checks out, 3 is a factor!
Because 3 is a factor of 18.
18 was unmarked while going through multiples when 3 was the number the algorithm was unmarking multiples for!
Start unmarking at the square
TL;DR: Start unmarking multiples of a number at that number squared.
All multiples below are already unmarked, because math is cool.
During the step the algorithm unmarks all multiples of a number.
We can start unmarking at that number squared.
Every smaller multiple was already unmarked in a previous iteration.
Why?
A multiple can be written as a multiplier times a number.
- m = multiple
- k = multiplier
- p = prime
m = k * p
The number that is now p, was previously k for every smaller prime number.
Because k * p = p * k, every multiple smaller than p * p has already been unmarked in a previous iteration.
Example
If our current detected prime, p = 5.
5 was previously the multiplier for every smaller prime number.
- 5 * 2 was unmarked when p was 2, we don't need to calculate 2 * 5
- 5 * 3 was unmarked when p was 3, we don't need to calculate 3 * 5
Step by step in code
The goal is to write a function that returns a list of prime numbers, up to upper_bound.
I named that variable
optimus_primewhile writing code for this post, because I thought that was funny
We initialise a list of booleans that is 1 bigger than the given upper_bound and call it sieve.
These booleans tell us if the number at that index is prime or not. (True for prime, False for not)
def primes_up_to(upper_bound):
# initialise sieve that marks all numbers as prime
sieve = [True] * (upper_bound + 1)
Smart people decided programmers start counting at 0, so that's why that list is 1 bigger than upper_bound.
It's also the reason why we have to unmark the index 0 along with the index 1 before we start our loop.
def primes_up_to(upper_bound):
# initialise sieve that marks all numbers as prime
sieve = [True] * (upper_bound + 1)
# 0 and 1 are not prime
sieve[0] = False
sieve[1] = False
This works out perfectly, because now every index exactly matches the number it represents.
You want to know if the number 69 is prime? The boolean at index 69 will tell you. Nice!
Loop over every number, starting at 2 and ending at the square root of upper_bound.
Inside the loop, index sieve with that number.
import math
def primes_up_to(upper_bound):
# initialise sieve that marks all numbers as prime
sieve = [True] * (upper_bound + 1)
# 0 and 1 are not prime
sieve[0] = False
sieve[1] = False
# iterate up to square root of upper_bound
# reason: if one factor of num is bigger than sqrt(upper_bound),
# an other factor _must_ be smaller than sqrt(upper_bound)
for num in range(2, math.floor(math.sqrt(upper_bound)) + 1):
# if sieve[num] is true, then num is prime
if sieve[num]:
If the boolean at that location is True, the number is prime and we unmark every multiple before moving on to the next step of our loop.
Do this by skip counting.
Start at the number squared and add the number until you hit upper_bound.
For every encountered multiple, set sieve at that number's index to False.
import math
def primes_up_to(upper_bound):
# initialise sieve that marks all numbers as prime
sieve = [True] * (upper_bound + 1)
# 0 and 1 are not prime
sieve[0] = False
sieve[1] = False
# iterate up to square root of upper_bound
# reason: if one factor of num is bigger than sqrt(upper_bound),
# an other factor _must_ be smaller than sqrt(upper_bound)
for num in range(2, math.floor(math.sqrt(upper_bound)) + 1):
# if sieve[num] is true, then num is prime
if sieve[num]:
# unmark all multiples
# start unmarking at num squared
# every smaller multiple has already been unmarked in previous iterations
for multiple in range(num ** 2, upper_bound + 1, num):
sieve[multiple] = False
At the end of the outer loop, sieve will be full of booleans corresponding to the primeness of every possible index to that list.
Use your favourite method to loop over a list while also getting the index, put the indexes with a true into a new list, and presto, primes.
import math
def primes_up_to(upper_bound):
# initialise sieve that marks all numbers as prime
sieve = [True] * (upper_bound + 1)
# 0 and 1 are not prime
sieve[0] = False
sieve[1] = False
# iterate up to square root of upper_bound
# reason: if one factor of num is bigger than sqrt(upper_bound),
# an other factor _must_ be smaller than sqrt(upper_bound)
for num in range(2, math.floor(math.sqrt(upper_bound)) + 1):
# if sieve[num] is true, then num is prime
if sieve[num]:
# unmark all multiples
# start unmarking at num squared
# every smaller multiple has already been unmarked in previous iterations
for multiple in range(num ** 2, upper_bound + 1, num):
sieve[multiple] = False
# sieve is done, turn `True` into numbers
return [idx for idx, mark in enumerate(sieve) if mark]
The returned value is a list of prime numbers, starting at 2, and ending in the last prime up to upper_bound.
-
primes_up_to(16)returns[2, 3, 5, 7, 11, 13]. -
primes_up_to(17)returns[2, 3, 5, 7, 11, 13, 17]. -
primes_up_to(18)returns[2, 3, 5, 7, 11, 13, 17]. -
primes_up_to(19)returns[2, 3, 5, 7, 11, 13, 17, 19].
Final code
pub fn primes_up_to(upper_bound: usize) -> Vec<usize> {
// initialise sieve that marks all numbers as prime
let mut sieve = vec![true; upper_bound + 1];
// 0 and 1 are not prime
sieve[0] = false;
sieve[1] = false;
// iterate up to square root of upper_bound
// reason: if one factor of num is bigger than sqrt(upper_bound),
// an other factor _must_ be smaller than sqrt(upper_bound)
for num in 2..=(upper_bound as f64).sqrt() as usize + 1 {
// if sieve[num] is true, then num is prime
if sieve[num] {
// unmark all multiples
// start unmarking at num squared
// every smaller multiple has already been unmarked in previous iterations
for multiple in (num * num..=upper_bound).step_by(num) {
sieve[multiple] = false;
}
}
}
// sieve is done, turn `true` into numbers
sieve
.iter()
.enumerate()
.filter_map(|(idx, mark)| match mark {
true => Some(idx),
false => None,
})
.collect()
}
function primesUpTo(upperBound) {
// initialise sieve that marks all numbers as prime
const sieve = Array.from({ length: upperBound + 1 }, () => true);
// 0 and 1 are not prime
sieve[0] = false;
sieve[1] = false;
// iterate up to square root of upperBound
// reason: if one factor of num is bigger than sqrt(upperBound),
// an other factor _must_ be smaller than sqrt(upperBound)
for (let num = 2; num <= Math.sqrt(upperBound) + 1; num++) {
// if sieve[num] is true, then num is prime
if (sieve[num]) {
// unmark all multiples
// start unmarking at num squared
// every smaller multiple has already been unmarked in previous iterations
for (let multiple = num ** 2; multiple <= upperBound; multiple += num) {
sieve[multiple] = false;
}
}
}
// sieve is done, turn `true` into numbers
const primes = [];
for (const [idx, mark] of sieve.entries()) {
mark && primes.push(idx);
}
return primes;
}
import math
def primes_up_to(upper_bound):
# initialise sieve that marks all numbers as prime
sieve = [True] * (upper_bound + 1)
# 0 and 1 are not prime
sieve[0] = False
sieve[1] = False
# iterate up to square root of upper_bound
# reason: if one factor of num is bigger than sqrt(upper_bound),
# an other factor _must_ be smaller than sqrt(upper_bound)
for num in range(2,math.floor(math.sqrt(upper_bound)) + 1):
# if sieve[num] is true, then num is prime
if sieve[num]:
# unmark all multiples
# start unmarking at num squared
# every smaller multiple has already been unmarked in previous iterations
for multiple in range(num**2, upper_bound + 1, num):
sieve[multiple] = False
# sieve is done, turn `True` into numbers
return [idx for idx, mark in enumerate(sieve) if mark]
Posted on August 13, 2021
Join Our Newsletter. No Spam, Only the good stuff.
Sign up to receive the latest update from our blog.