Sieve of Eratosthenes in Python

Why this matters

A prime number is a whole number greater than 1 whose only divisors are 1 and itself. The first primes are 2, 3, 5, 7, 11, 13 and so on. The number 4 is not prime, because 2 times 2 gives 4. A number that is not prime and is greater than 1 is called composite.

Many programs need a list of all the primes up to some limit, not just a test of one number. The simplest method is trial division: for each number, try every smaller divisor. That works, but it repeats a lot of work. The sieve of Eratosthenes is a method, named after the Greek mathematician Eratosthenes, that finds all the primes in one pass over a table. A sieve is a tool that lets small pieces through and holds back big ones. Here, the method holds back the composite numbers and lets the primes through.

This lesson builds the method one step at a time. It then writes the method in Python with a list of True and False values, counts how much work it does, and looks at the mistakes that come up in real code. A list is an ordered collection of items. An index is the position of an item in a list, and Python counts indexes from 0.

Big O notation describes how the work grows when the input grows. The lesson on time and space complexity in module 13 covers it in full. Here you only need the idea that the sieve does about n log log n steps for a limit n, which is much less work than trial division for large limits.

The idea with a small example

Take the limit 30. The goal is every prime below 30. Start with the numbers 2 to 29, none of them crossed out.

Step 1. Look at 2. Nothing before it has crossed it out, so 2 is prime. Cross out its multiples from 4 onward: 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26 and 28. Every even number is gone.

Step 2. Look at the next number that is not crossed out, which is 3. It is prime. Cross out its multiples from 9 onward: 9, 12, 15, 18, 21, 24 and 27. Some of these were crossed out by 2 already, and that is fine.

Step 3. The next number not crossed out is 5. It is prime. Its multiples from 25 onward are 25 and 30. The number 25 is new. The number 30 is not below the limit.

Step 4. The next number not crossed out is 7. It is prime. Its multiples from 49 onward are all at or above 30, so nothing is crossed out.

Step 5. The numbers 11, 13, 17, 19, 23 and 29 are all prime. Each one has a square at or above 30, so nothing new is crossed out.

The numbers that were never crossed out are 2, 3, 5, 7, 11, 13, 17, 19, 23 and 29. These are the primes below 30.

The same steps can be written as a Python program. It prints the numbers that get crossed out for each prime:

Save it as sieve_trace.py and run it with python3 sieve_trace.py. The output is:

p=2 marks: 4 6 8 10 12 14 16 18 20 22 24 26 28
p=3 marks: 9 12 15 18 21 24 27
p=5 marks: 25
p=7 marks: 
p=11 marks: 
p=13 marks: 
p=17 marks: 
p=19 marks: 
p=23 marks: 
p=29 marks: 
primes below 30: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

The line for 7 ends with the text marks: and nothing after it, because 7 times 7 is 49, which is past the limit.

Why the crossing out starts at p * p

The program starts crossing out at p * p, not at 2 * p. For p = 5, the multiples of 5 are 5, 10, 15, 20, 25 and so on. The numbers 10, 15 and 20 are already crossed out. Here is why. Any multiple of p that is smaller than p * p has the form k * p, where k is smaller than p. The number k has a prime factor that is smaller than p, and that smaller prime already crossed out k * p when it was processed. So the first multiple of p that might still be standing is p * p.

This saves writes. For p = 5, the program writes only 25 instead of 10, 15, 20 and 25. For large limits, the saving is large.

Why the outer loop can go all the way to the limit

The outer loop checks every number from 2 up to the limit. For a prime p with p * p at or above the limit, the inner loop does nothing. The outer loop still does one check per number, which costs very little. A version that stops the outer loop at the square root of the limit gives the same result. The version in this lesson keeps the loop simple. A shorter version that stops at the square root appears later in this lesson.

Writing it in Python

A list of True and False values is the sieve. Here is how it works.

[False] * n makes a list of n slots, with indexes 0 to n - 1. Every slot starts as False. The slot at index i means "i has been crossed out, so it is not prime" when it is True. The index is the number itself, so the list for a limit of 30 has 30 slots, and index 29 is the last one.

range gives the numbers that the inner loop visits. range(start, stop, step) gives start, start + step, and so on. It stops before stop, so stop is never included. For example, range(9, 30, 3) gives 9, 12, 15, 18, 21, 24 and 27.

Here is a complete program. It counts the primes below a limit and also counts the number of writes into the list, which is the amount of work done. The program is named sieve_count.py:

count_primes_below returns two values at once, written as a pair separated by a comma. Python calls this pair a tuple. The line count, marks = count_primes_below(n) takes the two values apart. The output is:

n        primes  marks     marks/(n*ln(ln n))
10       4       4         0.480
100      25      102       0.668
1000     168     1409      0.729
10000    1229    16979     0.765
100000   9592    193076    0.790
1000000  78498   2122046   0.808
 
Counting 0 and 1 as prime, n=10: 6
Correct count below 10: 4
 
Marks for n=1000000 starting at p*p: 2122046
Marks for n=1000000 starting at 2*p: 2775208

Read the code in three parts.

The outer loop runs p from 2 up to n - 1. If composite[p] is False, p was never crossed out, so p is prime. The program adds one to count.

The inner loop starts at p * p. range(p * p, n, p) visits p * p, then p * p + p, and so on, stopping before n. Each visit crosses out one multiple of p. Python integers never overflow, so p * p is always the exact product. The common mistakes section shows what that means.

The marks counter goes up once for every write. The table in the output compares the number of writes with n times the natural log of the natural log of n. The natural log, written ln, is the logarithm with base about 2.718. The next section explains what the comparison means.

Counting the work: n log log n

Each prime p writes into about n / p slots, because its multiples from p * p up to n are spaced p apart. The total work is the sum of n / p over all primes p up to n. That is n times the sum of 1 / p over the primes. The sum of 1 / p over the primes grows like log log n, which is the logarithm of the logarithm of n. The sum grows very slowly. For n = 1000000, the natural log of the natural log is about 2.6. So the number of writes is about n log log n.

The last column of the output shows the count of writes divided by n times ln(ln n). The ratio stays between 0.48 and 0.81 for all six limits. It rises slowly, and it does not grow in a way that would break the bound. This is evidence from six runs, not a proof, but it matches the idea that the work is a fixed fraction of n log log n. Unlike a timing, these counts are the same on every machine.

The outer loop does one check for each number from 2 to n - 1, which is n checks. That is also O(n). The total is O(n log log n), and for every practical limit the log log part is small.

The list holds one slot for each number below n. So the memory use is O(n). Trial division uses almost no memory, but it takes much more time for large limits. The sieve trades memory for speed.

The count column shows the number of primes below each limit. The values 4, 25, 168, 1229, 9592 and 78498 are the known counts of primes below 10, 100, 1000, 10000, 100000 and 1000000. Checking results against known values is a good habit when you write a new algorithm.

Second example: twin primes

Twin primes are two primes that differ by 2, such as 3 and 5, or 11 and 13. To find them, build the sieve once and then look at neighbouring pairs of numbers. The sieve is used as a lookup table. composite[i] answers the question in one step: False means i is prime.

This version also sets indexes 0 and 1 to True, because 0 and 1 are not prime. The table here includes the number 50, which is different from count_primes_below above. That function stops before n.

Output:

twin primes up to 50: (3, 5) (5, 7) (11, 13) (17, 19) (29, 31) (41, 43)

The loop range(2, limit - 1) stops before limit - 1. So the last value of i is limit - 2, and i + 2 never goes past the last slot of the table.

Third example: the largest gap between primes

The gap between two consecutive primes is the difference between them. The gap between 2 and 3 is 1, and between 7 and 11 it is 4. This example finds the largest gap below 1000. It uses the same sieve, then keeps track of the previous prime while it walks through the list.

Output:

largest gap between consecutive primes below 1000: 20 (after 887)

The value previous is None at the start. None is Python's word for "no value", and here it means no prime has been seen yet. The check previous is not None skips the first prime, because there is nothing before it to measure a gap from. The comparison uses >, not >=, so the first gap of the largest size is kept.

This example combines three ideas: the sieve builds a table, a loop reads the table in order, and a variable remembers the last good value. The same pattern appears when you read a list and compare each item with the one before it.

A shorter version with slices

Python can cross out a whole run of multiples in one statement. A slice is a part of a list. The expression list[start::step] picks every step-th item from start. Assigning a list to a slice replaces those items. The line in the function below crosses out all multiples of p from p * p on. The list on the right side has the same length as the slice, which len(range(p * p, n, p)) gives, so the assignment works.

Output:

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
78498

math.isqrt(n - 1) gives the whole part of the square root of n - 1. Any composite number below n has a prime factor at or below its square root, so the outer loop can stop there. The function expects n to be 2 or more. The expression [i for i in range(n) if is_prime[i]] is a list comprehension. It builds the list of primes in one line.

Common mistakes

Counting 0 and 1 as primes

The numbers 0 and 1 are not prime. The mistake happens when a program counts every slot that is still False, starting from index 0. The count_unmarked_from_zero function in sieve_count.py does exactly that. The expression sum(1 for flag in composite if not flag) adds 1 for each slot that is False. Its output for a limit of 10 is:

Counting 0 and 1 as prime, n=10: 6
Correct count below 10: 4

The six unmarked slots are 0, 1, 2, 3, 5 and 7. The correct answer is four: 2, 3, 5 and 7. The fix is to start the count at 2, as the outer loop does, or to set composite[0] and composite[1] to True before the loop, as the twin primes example does.

Expecting an integer overflow

In Java, an int holds numbers only up to 2147483647, and a larger product wraps around to a wrong value. Python does not work that way. An int in Python grows to hold any size, so p * p is always the exact product. The program below shows this:

Output:

4294967296
1267650600228229401496703205376
2147488281

The expression 2 ** 100 means 2 raised to the power 100. Because Python does not overflow, Java would need a long cast for this product, and Python has no such cast.

Index out of range

A list of length 30 has indexes 0 to 29. Writing to index 30 raises an error. This program shows the message Python prints:

composite = [False] * 30
composite[30] = True

Running it prints a traceback. The traceback below is shown with the folder names replaced by .... The file name and the line number are the ones Python prints:

Traceback (most recent call last):
  File ".../sieve_index_error.py", line 2, in <module>
    composite[30] = True
    ~~~~~~~~~^^^^
IndexError: list assignment index out of range

The line with ^ marks the part of the code that failed, which is composite[30]. The message says the assignment failed, but it does not show the index or the length. The fix is to make sure the loop stops at m < n, not m <= n, when the list has n slots. If you want primes up to and including n, make the list n + 1 slots long, as the twin primes example does.

Negative indexes do not raise an error

Python accepts negative indexes. composite[-1] means the last slot, which is index 29 for a list of 30 slots. A bug that makes an index negative writes to the wrong slot, and Python does not print an error. The program below shows this:

Output:

True
True

If you know Java, this is different from Java, which raises an error for an index of -1.

Using < where you need <=

count_primes_below(7)[0] is 3, because it counts 2, 3 and 5. It does not count 7, because 7 is not below 7. count_primes_below(8)[0] is 4. The difference is the word "below". Read the task carefully to see whether the limit is included.

Starting the inner loop at 2 * p

Starting the inner loop at 2 * p gives the right answer, because the numbers between 2 * p and p * p are multiples of p that were already crossed out. The program still does more writes. The output above shows 2775208 writes for 2 * p and 2122046 for p * p, both for a limit of 1000000. The results are the same, but the first version does about 31 percent more work.

A step of zero

If the outer loop starts at p = 0, the inner loop begins at 0 * 0, which is 0, and the step is 0. A range with a step of 0 is not allowed. The program below shows what Python does:

composite = [False] * 30
p = 0
for m in range(p * p, len(composite), p):
    composite[m] = True

Running it prints a traceback, shown here with the folder names replaced by ...:

Traceback (most recent call last):
  File ".../sieve_zero_step.py", line 3, in <module>
    for m in range(p * p, len(composite), p):
             ~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: range() arg 3 must not be zero

A ValueError is an error that Python raises when a function gets a value of the right type with a wrong value. If you know Java, this is different: in Java the same mistake makes the loop run forever. The cause is the step: any prime used in the step must be at least 2.

Practice

Exercise 1: count the writes for a limit of 20

Using the method from this lesson, work out how many writes the inner loop makes for n = 20. Count writes for each prime p, not the number of distinct numbers crossed out.

Show the answer

For p = 2, the writes are at 4, 6, 8, 10, 12, 14, 16 and 18. That is 8 writes. For p = 3, the writes are at 9, 12, 15 and 18. That is 4 writes. For p = 5, the first multiple is 25, which is not below 20, so there are no writes. The total is 8 + 4 = 12. Running count_primes_below(20)[1] from sieve_count.py gives the same 12.

Exercise 2: list the primes up to and including 50

Write a function that returns the primes from 2 up to and including 50, and print the list.

Show the answer

There are 15 primes: 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43 and 47. The number 50 is not prime, so including it does not change the list. Call primes_below(51) from the slices example, which returns every prime below 51. Printing it gives [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47].

Exercise 3: find the twin primes below 50

Use the table from the twin primes example to list every pair of primes that differ by 2, with both numbers below 50.

Show the answer

The pairs are (3, 5), (5, 7), (11, 13), (17, 19), (29, 31) and (41, 43). The number 47 has no partner, because 45 and 49 are not prime.

Exercise 4: explain why the crossing out starts at p * p for p = 5

Explain in plain words why, for p = 5, the crossing out can start at 25 and skip 10, 15 and 20.

Show the answer

The number 10 is 2 times 5, so 2 crossed it out when 2 was processed. The number 15 is 3 times 5, so 3 crossed it out. The number 20 is 2 times 10, so 2 crossed it out too. The first multiple of 5 that no smaller prime has touched is 25, which is 5 times 5. The same reasoning applies to every prime p: any smaller multiple of p has a smaller prime factor.

Exercise 5: predict the product in Python

Without running the code, predict the output of this program. Then run it to check.

p = 46341
print(p * p)

Show the answer

The output is 2147488281. In Java, the same product would wrap around to -2147479015, because it is larger than the largest int. Python does not wrap around, so the output is the exact product.

Summary

• A prime has no divisors other than 1 and itself. The sieve finds all primes below a limit in one pass.
• For each prime p, cross out its multiples starting at p * p, because smaller multiples were crossed out earlier.
• A list of True and False values, with one slot per number, is the table, and False means not yet crossed out.
• The work is about n log log n writes, which is far less than checking each number by trial division.
• Check the list indexes, remember that range stops before its end, and do not count 0 and 1 as prime.