The versatile sieve of Eratosthenes
Introduction
The sieve of Eratosthenes is an efficient algorithm to mark whether each integer between 0 to n is prime. It turns out that the algorithm is useful for calculating other interesting values that are also based on the prime factorization of a number.
This page will show how the basic sieve of Eratosthenes can be modified to calculate other number-theoretic functions. A prior understanding of the basic algorithm will be assumed (please browse other resources for an introduction).
Designing and applying these modified forms of the sieve of Eratosthenes has been very helpful in a number of my Project Euler solution implementations.
Technical notes
The sieve of Eratosthenes (and its variants) takes an upper limit integer n, and computes a table of some particular number-theoretic function f(k) for every value of k from 0 to n. Hence if you wish to evaluate f(k) for many values of k between 0 and n (i.e. the usage is dense), then it helps to run a sieve first to pre-compute all the values in that domain. If the queries of k are sparse however (e.g. only testing if a couple of numbers are prime), then sieving is not efficient in time or memory.
All code examples are given in the Python programming language, for clarity and to avoid issues with integer overflow. These algorithms can certainly be implemented in languages with fixed-width integer types (such as C/C++/C#/Java/JavaScript/Rust), but extreme care must be taken to prove that every calculation does not overflow the chosen integer type’s range.
Each example sieving function takes a non-negative integer named limit
and returns a list/array of limit+1
elements (usually integers). Then, for any query number k
in the range 0 <= k <= limit
, the expression result[k]
gives the value of the desired number-theoretic function.
Most number-theoretic functions don’t have a meaningful or consistent definition at 0 or 1. But because Python arrays start at index 0, I made an earnest attempt to give sensible values to indexes 0 and 1. For the most part, you can ignore these first two values, and only consider values starting at index 2.
Source code
- Python (plain)
- Python (NumPy)
- Java
-
- EratosthenesSieves.java (SE 1.2+)
- EratosthenesSievesTest.java (SE 8+)
- TypeScript / JavaScript
- C#
- C++ (C++11 and above)
- C (C99 and above)
- Rust
License: Project Nayuki hereby places all code on this page regarding the sieve of Eratosthenes in the public domain. Retaining the credit notice containing the author and URL is encouraged but not required.
Primeness (the basic sieve)
def sieve_primeness(limit): result = [True] * (limit + 1) result[0] = False result[1] = False for i in range(2, len(result)): if result[i]: for j in range(i * i, len(result), i): result[j] = False return result
isPrime(k) is true if k is a prime number; otherwise it’s false if k is 0 or 1 or composite.
Table of example values:
k | isPrime(k) |
---|---|
0 | false |
1 | false |
2 | true |
3 | true |
4 | false |
5 | true |
6 | false |
7 | true |
8 | false |
9 | false |
10 | false |
11 | true |
12 | false |
13 | true |
14 | false |
15 | false |
16 | false |
17 | true |
18 | false |
19 | true |
20 | false |
21 | false |
22 | false |
23 | true |
24 | false |
25 | false |
26 | false |
27 | false |
28 | false |
29 | true |
30 | false |
Smallest prime factor
def sieve_smallest_prime_factor(limit): result = [0] * (limit + 1) result[1] = 1 for i in range(2, len(result)): if result[i] == 0: result[i] = i for j in range(i * i, len(result), i): if result[j] == 0: result[j] = i return result
smallestPrimeFactor(k) is the smallest integer p such that p is greater than 1 and p divides k. The smallest proper divisor is necessarily a prime number.
When given a composite number k, we can repeatedly divide out its smallest prime factor in order to obtain its prime factorization.
Table of example values:
k | smallestPrimeFactor(k) |
---|---|
0 | 0 |
1 | 1 |
2 | 2 |
3 | 3 |
4 | 2 |
5 | 5 |
6 | 2 |
7 | 7 |
8 | 2 |
9 | 3 |
10 | 2 |
11 | 11 |
12 | 2 |
13 | 13 |
14 | 2 |
15 | 3 |
16 | 2 |
17 | 17 |
18 | 2 |
19 | 19 |
20 | 2 |
21 | 3 |
22 | 2 |
23 | 23 |
24 | 2 |
25 | 5 |
26 | 2 |
27 | 3 |
28 | 2 |
29 | 29 |
30 | 2 |
Notice that for k > 1, smallestPrimeFactor(k) = k if and only if k is a prime number. Hence, this information can be viewed as a richer version of the plain sieve of Eratosthenes (which only gives a Boolean value per number).
Totient function
def sieve_totient(limit): result = list(range(limit + 1)) for i in range(2, len(result)): if result[i] == i: for j in range(i, len(result), i): result[j] -= result[j] // i return result
totient(k) (or φ(k)) is defined as the number of integers in the range [1, k] that are coprime with k. This value is useful for modular arithmetic.
Table of example values:
k | totient(k) |
---|---|
0 | 0 |
1 | 1 |
2 | 1 |
3 | 2 |
4 | 2 |
5 | 4 |
6 | 2 |
7 | 6 |
8 | 4 |
9 | 6 |
10 | 4 |
11 | 10 |
12 | 4 |
13 | 12 |
14 | 6 |
15 | 8 |
16 | 8 |
17 | 16 |
18 | 6 |
19 | 18 |
20 | 8 |
21 | 12 |
22 | 10 |
23 | 22 |
24 | 8 |
25 | 20 |
26 | 12 |
27 | 18 |
28 | 12 |
29 | 28 |
30 | 8 |
Notice that totient(k) = k−1 if and only if k is a prime number. Hence, this information can be viewed as a richer version of the plain sieve of Eratosthenes (which only gives a Boolean value per number).
Omega function
def sieve_omega(limit): result = [0] * (limit + 1) for i in range(2, len(result)): if result[i] == 0: for j in range(i, len(result), i): result[j] += 1 return result
omega(k) (or ω(k)) is defined as the number of distinct prime factors that k has.
Table of example values:
k | omega(k) |
---|---|
0 | 0 |
1 | 0 |
2 | 1 |
3 | 1 |
4 | 1 |
5 | 1 |
6 | 2 |
7 | 1 |
8 | 1 |
9 | 1 |
10 | 2 |
11 | 1 |
12 | 2 |
13 | 1 |
14 | 2 |
15 | 2 |
16 | 1 |
17 | 1 |
18 | 2 |
19 | 1 |
20 | 2 |
21 | 2 |
22 | 2 |
23 | 1 |
24 | 2 |
25 | 1 |
26 | 2 |
27 | 1 |
28 | 2 |
29 | 1 |
30 | 3 |
Note that for a given m, the smallest k such that omega(k) = m is given by k = primorial(m). Furthermore the function grows rather slowly, because omega(k) ≤ log2(k).
Radical function
def sieve_radical(limit): result = [1] * (limit + 1) result[0] = 0 for i in range(2, len(result)): if result[i] == 1: for j in range(i, len(result), i): result[j] *= i return result
radical(k) is defined as the product of the unique prime factors of k.
Table of example values:
k | radical(k) |
---|---|
0 | 0 |
1 | 1 |
2 | 2 |
3 | 3 |
4 | 2 |
5 | 5 |
6 | 6 |
7 | 7 |
8 | 2 |
9 | 3 |
10 | 10 |
11 | 11 |
12 | 6 |
13 | 13 |
14 | 14 |
15 | 15 |
16 | 2 |
17 | 17 |
18 | 6 |
19 | 19 |
20 | 10 |
21 | 21 |
22 | 22 |
23 | 23 |
24 | 6 |
25 | 5 |
26 | 26 |
27 | 3 |
28 | 14 |
29 | 29 |
30 | 30 |
Notice that radical(k) = k if and only if k is square-free.