## DEV Community is a community of 752,525 amazing developers

We're a place where coders share, stay up-to-date and grow their careers.

Nicky Meuleman

Posted on • Originally published at nickymeuleman.netlify.app

# The Sieve of Eratosthenes

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.

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_prime` while 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 = False
sieve = 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 = False
sieve = 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 = False
sieve = 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 = False
sieve = 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 = false;
sieve = 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 = false;
sieve = 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 = False
sieve = 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]
``````

## Discussion (4) chenge
``````pub fn primes_up_to(limit: u32) -> Vec<u32> {
let mut primes = Vec::new();
let mut candidates: Vec<_> = (2..=limit).rev().collect();
while let Some(prime) = candidates.pop() {
primes.push(prime);
candidates.retain(|n| n % prime != 0);
}
primes
}
``````