Sieve of Eratosthenes

Why this matters

A prime number is a whole number greater than 1 that has no divisors except 1 and itself. The numbers 2, 3, 5, 7, 11 and 13 are prime. The numbers 4, 6, 8, 9 and 10 are not. Each of them can be written as a product of smaller whole numbers: 4 is 2 times 2, 6 is 2 times 3, and so on. A number like that is called composite.

Primes appear in encryption and in many math problems that split numbers into parts. Interview tasks often ask for every prime up to a limit, such as 1,000 or 1,000,000. You can answer that question in more than one way. This lesson teaches the fastest common method for a list of primes, the Sieve of Eratosthenes. The Greek mathematician Eratosthenes described it about 2,200 years ago. A sieve is a kitchen strainer: the things you want stay in, and the things you do not want fall through. Here the unwanted numbers are the composites, and the program crosses them out.

The code uses for loops, if, and continue, which you met earlier in this module. It also uses std::vector<bool>, a list of true and false values. The std:: prefix marks a name from the standard library. The program in this lesson is written for C++20. You can compile any sample with g++ -std=c++20 -Wall.

The slow way: test each number

Start with the method you would try first. To check whether 13 is prime, try each whole number from 2 up to 12 and see whether any of them divides 13 exactly. A divisor of a number is a number that divides it with no remainder. In C++ the % operator gives the remainder after division, so 13 % 2 is 1, and 12 % 2 is 0. If m % d is 0, then d is a divisor of m.

For 13, none of the numbers from 2 to 12 divides it, so 13 is prime. For 12, the number 2 divides it, so the check stops after one try. This method is called trial division. It is correct, but it does a lot of work when the limit is large. For a number m, it may need up to m - 2 checks, and you repeat that for every number up to the limit.

The sample in the sieve-counting section below measures this. For limit 10,000 the simple method makes about 5.8 million divisor checks. The sieve makes about 17,000 writes for the same limit. The gap grows as the limit grows.

The idea of the sieve

The sieve does not test numbers one at a time. It starts by assuming every number is prime. Then it crosses out the numbers that it can prove are not prime. A composite number always has a prime divisor that is smaller than it. So when the program reaches a prime p, it crosses out every multiple of p that is still in the list.

Here are the steps for the limit 30. Each step is one thing the program does.

1. Write down the numbers 2 to 30. Every number is a candidate, which means it might still be prime.
2. Look at the first candidate, 2. It has not been crossed out, so it is prime. Cross out its multiples: 4, 6, 8, 10, and so on up to 30.
3. Look at the next candidate, 3. It has not been crossed out, so it is prime. Cross out its multiples: 6, 9, 12, 15, and so on up to 30. Some of these were already crossed out in step 2, and crossing them out again changes nothing.
4. The next number is 4. It was crossed out in step 2, so it is skipped.
5. The next candidate is 5. It is prime. Cross out its multiples: 10, 15, 20, 25 and 30. Most were already crossed out.
6. The next candidate is 7. Its square, 49, is larger than 30, so the sieve stops. The next part explains why that is safe.
7. Every number still not crossed out is prime.

Worked example 1: a trace for the limit 30

Here is the trace in the form the program follows. A trace is a record of what happens at each step, written down by hand.

• Start with the numbers 2 to 30 as candidates.
• p = 2 is a candidate, so 2 is prime. Mark 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28 and 30. That is 14 marks.
• p = 3 is a candidate, so 3 is prime. Mark 9, 12, 15, 18, 21, 24, 27 and 30. That is 8 marks. Four of them, 12, 18, 24 and 30, were already marked.
• p = 4 is not a candidate, so skip it.
• p = 5 is a candidate, so 5 is prime. Mark 25 and 30. That is 2 marks. The multiples 10, 15 and 20 were already marked.
• p = 6 is not a candidate. The program has also reached the end of the useful range, because 6 times 6 is 36, which is larger than 30.

The unmarked numbers are 2, 3, 5, 7, 11, 13, 17, 19, 23 and 29. These are the ten primes up to 30. The total number of marks written was 14 + 8 + 2 = 24, even though only 19 numbers from 2 to 30 are composite. Some composites were marked more than once, such as 30, which was marked by 2, 3 and 5.

Why start at p times p

Look at the step for p = 5 again. The multiples 10, 15 and 20 were already crossed out. The first multiple of 5 that is still a candidate is 25, which is 5 times 5.

The same reason holds for every prime. When the program reaches a prime p, each smaller multiple of p equals p times some number k, where k is smaller than p. The number k has a prime factor that is smaller than p, and that factor already crossed out the multiple. So the program can begin crossing out at p times p.

Now the stopping rule. Every composite number m has a prime divisor p with p times p less than or equal to m. For example, 30 has the prime divisor 2, and 2 times 2 is 4, which is at most 30. So once p times p is larger than the limit, every composite number up to the limit has already been crossed out by a smaller prime, and the loop can end.

With a limit of 30, the primes that need to do any work are 2, 3 and 5, because 5 times 5 is 25 and 7 times 7 is 49. The loop stops once the square of p is past the limit.

The first program

The program below sieves the numbers up to 30 and prints the primes. The code is written so that the index of the vector is the number itself. Index 7 holds the answer for the number 7. A vector is a resizable list from the standard library. An index is a position in the list, counted from 0.

Output:

2 3 5 7 11 13 17 19 23 29

Read the program in four parts.

• std::vector<bool> is_prime(n + 1, true); creates n + 1 boxes, numbered 0 to n. Each box holds true at the start, which means "still a candidate". The n + 1 is needed because index n must exist.
• The two lines with = false clear 0 and 1, which are not prime.
• The outer loop runs p from 2 while p * p <= n. The inner loop starts at p * p and adds p each time. m += p means "add p to m". So the inner loop visits p times p, then p times p plus p, and so on, until m passes n.
• The last loop reads the boxes from 2 to n and prints each index whose box still says true.

The continue in the outer loop skips the rest of the loop body and moves on to the next value of p. It is used when p was already crossed out, so p is not prime and needs no work.

Worked example 2: count the primes below a limit

Often you need the number of primes, not the list. The next program counts the primes that are smaller than n. The word "below" matters here. Below 10 means smaller than 10, so 10 itself is not counted.

Output:

primes below 10: 4
primes below 30: 10
primes below 100: 25
primes below 1000: 168
primes below 1000000: 78498

Three details are new here.

First, the vector has exactly n boxes, indexes 0 to n - 1, so every box is for a number smaller than n. The loops use < instead of <= for the same reason. Using < for every loop bound keeps the code consistent. The idea of a range that includes its start but stops before its end has its own lesson in this course.

Second, the guard if (n <= 2) return 0; handles small limits. For n equal to 1, the vector has one box, index 0, and the line is_prime[1] = false; writes outside the vector. The "Common mistakes" section below shows what that looks like.

Third, the counting loop is separate from the sieve. The sieve marks the boxes, and the second loop counts the boxes that still say true. Keeping the two jobs apart makes each loop easier to check.

Worked example 3: smallest prime factor

The sieve can store more than yes or no. If the program stores the smallest prime that divides each number, it can split any number into primes quickly. This is useful when a program must factor many numbers, for example in a number-theory exercise.

Output:

84 = 2 2 3 7
97 = 97
2 = 2
1 = (no prime factors)

The table smallest has one entry for each number. An entry is 0 while the number is still a candidate. The first prime to reach a multiple writes its own value into the entry, and later primes leave it alone, so each entry holds the smallest prime factor. The loop if (smallest[p] != 0) continue; skips composite numbers as before, because a composite was already written by a smaller prime.

To factor 84, the program reads smallest[84], which is 2. It prints 2 and divides 84 by 2 to get 42. Then it reads smallest[42], which is 2, and so on. The result is 2, 2, 3 and 7, because 2 times 2 times 3 times 7 is 84.

This version uses one int per number. For a limit of ten million, that is about 40 MB, which is the number of bytes in ten million 4-byte values. A vector<bool> version uses one bit per number instead, but it cannot store the prime factor. Choose the table that holds the information you need.

How much work the sieve does

A program's work is often measured by counting the basic steps it performs. Two steps are counted here. One is a write, which means setting a box to false in the sieve. The other is a divisor check, which means one % test in trial division. The next program counts both for four limits.

Output:

n=30  sieve writes=24  trial division checks=135
n=100  sieve writes=104  trial division checks=1133
n=1000  sieve writes=1411  trial division checks=78022
n=10000  sieve writes=16981  trial division checks=5775223

The trial division count grows very fast. From 1,000 to 10,000, the limit grows ten times, and the checks grow about 74 times. The sieve writes grow about 12 times in the same range.

The sieve's count is in the style of O(n log log n). In O notation, the letter O describes how the work grows as the input, here n, grows. Log log n means the logarithm of the logarithm of n. For a prime p, the sieve writes about n divided by p times. Adding these counts over the primes up to the square root of n gives about n times a sum of fractions 1/2 + 1/3 + 1/5 + 1/7 and so on, and that sum grows very slowly. For any limit that fits in memory, log log n is a small number, so the sieve does a little more than n writes. It is close to linear.

The count of writes can be larger than the count of composites, as the trace showed. The sieve writes 24 times for the limit 30, but only 19 numbers from 2 to 30 are composite. The extra writes are numbers that were crossed out more than once. This does not make the sieve wrong. It only means the marking work is not exactly one write per composite.

Common mistakes

Each mistake below is one that a beginner makes often. Each output shown came from a test run of code that has the mistake in it.

Forgetting that 0 and 1 are not prime

If the program does not clear the boxes for 0 and 1, they stay true, and they are counted as primes.

std::vector<bool> is_prime(n, true);  // forgot: is_prime[0] = false; is_prime[1] = false;

Compiling a version of count_primes_below with that line removed and calling it with 10 prints:

6

The correct answer is 4. The two extra counts are 0 and 1. The fix is the two lines that set is_prime[0] and is_prime[1] to false, as the code above does.

Writing past the end of the vector for small limits

A std::vector does not check indexes by default. If you write to an index that does not exist, the program may crash, print a wrong answer, or seem to work. The line is_prime[1] = false; is a problem when the limit is 1, because the vector then has only index 0.

Here is a complete program with that bug:

#include <vector>
 
int main() {
    int n = 1;
    std::vector<bool> is_prime(n, true);
    is_prime[0] = false;
    is_prime[1] = false;  // index 1 is past the end when n is 1
}

Compile it with the C++ library's own checks turned on, so the bad write is reported:

g++ -std=c++20 -Wall -D_GLIBCXX_ASSERTIONS -g -o bad2 bad2.cpp
./bad2

The program prints this message (the path and line number depend on your GCC version) and stops with exit status 134, the code the system reports when a program stops itself after a failed check:

/usr/include/c++/15/bits/stl_bvector.h:1158: constexpr std::vector<bool, _Alloc>::reference std::vector<bool, _Alloc>::operator[](size_type) [with _Alloc = std::allocator<bool>; reference = std::vector<bool>::reference; size_type = long unsigned int]: Assertion '__n < this->size()' failed.

The message says the index was not less than the size. The fix is the guard if (n <= 2) return 0;. Without -D_GLIBCXX_ASSERTIONS, the same bug may pass silently, so turn the checks on while you test. The flag works with GCC and with the libstdc++ library that ships with it.

A loop bound that is one too large

The loop below uses <= with a vector of size n. Index n does not exist.

std::vector<bool> is_prime(n, true);   // valid indexes are 0 to n - 1
for (int m = 0; m <= n; ++m) {         // reaches index n
    is_prime[m] = false;
}

With -D_GLIBCXX_ASSERTIONS and n equal to 10, the program stops with the same kind of message as above, with the same Assertion '__n < this->size()' failed. text. The fix is to use < for a vector of size n, or to make the vector n + 1 and keep <=. Choose one style and use it in every loop of the program.

Overflow in the square of p

The loop test p * p <= n multiplies two int values. Overflow means a result is too big for the space the type gives it. In C++, signed integer overflow is undefined behaviour, which means the program is no longer guaranteed to do anything sensible. The example below uses 50,000, which is small, to show the problem.

#include <iostream>
 
int main() {
    int n = 50000;
    int square = n * n;   // 2500000000 does not fit in int
    std::cout << square << '\n';
}

Compile with the undefined-behaviour checker turned on:

g++ -std=c++20 -Wall -fsanitize=undefined -g -o bad4 bad4.cpp
./bad4

The output is:

bad4.cpp:5:9: runtime error: signed integer overflow: 50000 * 50000 cannot be represented in type 'int'
-1794967296

A common int holds values up to 2147483647. For a sieve with a limit near two billion, the square of p can pass that size. Store p and the limit in long long in that case. For the limits in this lesson, int is enough.

Skipping the check for crossed-out numbers

The continue line in the outer loop is not needed for a correct answer. Without it, the program also marks multiples of composite numbers. The primes are still right, because each composite multiple was already marked by a prime. The program simply does more writes.

for (int p = 2; p * p <= n; ++p) {
    // no "if (!is_prime[p]) continue;" line
    for (int m = p * p; m <= n; m += p) {
        is_prime[m] = false;
    }
}

Counting the writes for the same limits gives these numbers: 146 writes instead of 104 for n equal to 100, and 2550 instead of 1411 for n equal to 1000. The skip is cheap, and it keeps the work down.

Practice

Exercise 1: trace by hand

Trace the sieve for the limit 20. Write down every prime from 2 to 20, and count the writes made for each prime p.

Show the answer

The primes are 2, 3, 5, 7, 11, 13, 17 and 19. For p = 2, the program marks 4, 6, 8, 10, 12, 14, 16, 18 and 20, which is 9 writes. For p = 3, it marks 9, 12, 15 and 18, which is 4 writes. For p = 4, the program skips, because 4 is already marked. The loop stops before p = 5, because 5 times 5 is 25, which is larger than 20. The total is 9 + 4 = 13 writes. A test run of the same loops for the limit 20 also counts 13 writes.

Exercise 2: the largest prime below a limit

Write a function largest_prime_below(int n) that returns the largest prime smaller than n. Return 0 when no prime exists. Test it with 10, 100, 3 and 2.

Show the answer

int largest_prime_below(int n) {
    if (n <= 2) return 0;
    std::vector<bool> is_prime(n, true);
    is_prime[0] = is_prime[1] = false;
    for (int p = 2; p * p < n; ++p) {
        if (!is_prime[p]) continue;
        for (int m = p * p; m < n; m += p) is_prime[m] = false;
    }
    for (int i = n - 1; i >= 2; --i) {
        if (is_prime[i]) return i;
    }
    return 0;
}

The results are 7 for 10, 97 for 100, 2 for 3, and 0 for 2. The last loop runs down from n - 1, so the first prime it meets is the largest one.

Exercise 3: primes in a range

Write a program that prints every prime from 100 to 150. Use one sieve for the whole job.

Show the answer

Sieve up to 150, then print only the numbers from 100 to 150 that are still marked. The output is:

101 103 107 109 113 127 131 137 139 149 

The sieve does not need a second pass for each range, because the crossed-out boxes already answer the question for every number up to the limit.

Exercise 4: factor a number

Using the smallest-prime-factor table from worked example 3, print the prime factors of 360, with each factor printed once for each time it divides the number.

Show the answer

Build the table up to 360, then divide repeatedly. The output is:

2 2 2 3 3 5 

Check it: 2 times 2 times 2 times 3 times 3 times 5 is 360.

Exercise 5: predict the work

Without running the program, predict how many writes the sieve makes for the limit 10. Use the loop bounds from the first program, with <= and p * p <= n. Then run it to check.

Show the answer

For p = 2, the program marks 4, 6, 8 and 10, which is 4 writes. For p = 3, it marks 9, which is 1 write, because 3 times 3 is 9 and 12 is past the limit. The loop then stops, because 4 times 4 is 16, which is larger than 10, and 4 is not prime anyway. The total is 4 + 1 = 5 writes.

Summary

• The Sieve of Eratosthenes starts with every number as a candidate and crosses out the multiples of each prime.
• When the sieve reaches a prime p, it begins at p times p, and it stops once p times p passes the limit.
• The count of primes below n is the number of boxes that still say true after the sieve, for indexes from 2 to n - 1.
• The sieve does about n log log n writes, while trial division does far more divisor checks for the same limit.
• Store the smallest prime factor instead of true or false when you need to split numbers into primes.
• Check small limits, index bounds and int overflow when you test. The compiler flags -D_GLIBCXX_ASSERTIONS and -fsanitize=undefined report these bugs when they happen.