Why is my Sieve of Eratosthenes using generators so slow

Question

This is my prime sieve based on the Sieve of Eratosthenes.

import itertools
sieve = itertools.count(2)
primes = []
nprime = 0
for _ in range(10000):
    nprime = next(sieve)
    primes.append(nprime)
    sieve = filter(eval(f"lambda x: x % {nprime} != 0"), sieve)

On my computer this code takes 13.4 seconds to generate 10,000 primes.

Answer 1

That's not the sieve of Eratosthenes

It's rather prime trial division.

Your filter for 2 lets 1/2 of all numbers through, your filter for 3 lets 2/3 of all remaining numbers through, etc. So 1/2 * 2/3 * 4/5 * 6/7 = 22% of all numbers make it through your filters for 2, 3, 5 and 7, so your filter for 11 will have to check 22% of all numbers. The real sieve of Eratosthenes treats 11 by marking off every 11th number, i.e., only 9% of all numbers. Less than half of your percentage. And yours gets worse for the larger primes. For example your filters for the primes from 2 to 89 let ~12% of all numbers through, so your filter for 97 checks ~12% of all numbers. While the real sieve of Eratosthenes treats 97 by marking off only about 1% of all numbers.

So that's why it's slow. It's doing way too much work.

An improvement

You can make it a bit faster by not stacking filters but simply doing the trial divisions at the current number, using the primes you're collecting anyway:

from itertools import count

def primes(n):
    primes = []
    sieve = (x for x in count(2) if all(x % p for p in primes))
    for _ in range(n):
        primes.append(next(sieve))
    return primes

Or using islice:

from itertools import count, islice

def primes(n):
    primes = []
    sieve = (x for x in count(2) if all(x % p for p in primes))
    primes += islice(sieve, n)
    return primes

Demo:

>>> primes(25)
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

The real thing

But here's a genuine lazy one based on the paper mentioned by trentcl. I add each prime to the list of its next multiple. When the list of the next x is empty, it's thus a prime. Otherwise, I move the primes from the list to their next multiples.

from itertools import count
from collections import defaultdict

def primes(n):
    marked = defaultdict(list)
    primes = []
    for x in count(2):
        if x in marked:
            for p in marked.pop(x):
                marked[x + p].append(p)
        else:
            primes.append(x)
            if len(primes) == n:
                return primes
            marked[x * x].append(x)

Benchmark 1

Times for n = 10000 (benchmark code at the end of the answer):

9.330 s  original
3.502 s  improved1
3.499 s  improved2
0.065 s  genuine

That said, your original and my improved version can be made much faster by using fewer filters (only up to the square root of the current number), as I've done in my later answer.

To infinity and beyond

As Eric Duminil showed, it's trivial to make the latter an infinite generator, which also makes it shorter. No need to build a result list and check whether you're done, just yield the primes:

infinite:

def primes():
    marked = defaultdict(list)
    for x in count(2):
        if x in marked:
            for p in marked.pop(x):
                marked[x + p].append(p)
        else:
            yield x
            marked[x * x].append(x)

The other solutions can similarly be turned into infinite generators, but meh, they're no good anyway. Plus in my improved1 and improved2, the primes list is not just the result but is also used for filtering, so I need to build it anyway.

As infinite generator, it's more versatile. You can still get a list of the first n if you want to:

>>> list(islice(primes(), 20))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]

You can also get all primes up to a desired limit:

>>> *itertools.takewhile(72 .__gt__, primes()),
(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71)

Or if you neither know how many you need nor how big you need but you need "enough" somehow, for example enough so that their product exceeds a million:

>>> result, product = [], 1
>>> it = primes()
>>> while product < 10**6:
        p = next(it)
        result.append(p)
        product *= p

>>> result
[2, 3, 5, 7, 11, 13, 17, 19]

Or if you just want the millionth prime:

>>> next(islice(primes(), 999999, None))
15485863

Or if you have a prime and want to know the "how manyth" prime it is:

>>> from operator import indexOf
>>> indexOf(primes(), 15485863) + 1
1000000

Or if you have a prime and want to know the next larger one (granted, there are more efficient ways to do that):

>>> it = primes()
>>> 89 in it and next(it)
97

Note that in the last few cases, you don't want a list of primes at all, so building one is wasted memory.

Speaking of memory, though... the generator (and also my version) do take a lot of memory. Here are memory measurements and benchmarks for up to n=1,000,000 for the above "infinite" solution and the solutions I'll show below:

            |     time for the first...    | tracemalloc for 1,000,000
            |  10,000   100,000  1,000,000 |     size         peak
------------+------------------------------+--------------------------
   infinite | 0.068 s   0.971 s   12.450 s |  231,866,275  312,793,599
   no_lists | 0.035 s   0.574 s    8.065 s |  143,870,959  226,554,527
     lazier | 0.033 s   0.456 s    6.016 s |   36,500,174   36,535,742
  recursive | 0.032 s   0.473 s    6.116 s |       72,169      109,109
    twostep | 0.015 s   0.217 s    2.897 s |       71,986      108,926

The memory usage was measured like this:

import tracemalloc
tracemalloc.start()
it = primes()
prime = next(islice(it, 9999999, None))
size, peak = tracemalloc.get_traced_memory()

So the above "infinite" solution takes hundreds of megabytes to generate the first million primes, i.e., hundreds of bytes per prime. Quite a lot. Let's see why, by doing global marked first thing in the primes() function so we can inspect it afterwards. We'll see that almost a million composites are marked and the dictionary takes 84 MB:

>>> it = primes()
>>> next(islice(it, 999999, None))
15485863
>>> from sys import getsizeof
>>> len(marked), getsizeof(marked)
(999919, 83886176)

The lists of primes total 88 MB, and almost all of them contain only a single prime marking the composite:

>>> sum(map(getsizeof, marked.values()))
87992872
>>> from collections import Counter
>>> Counter(map(len, marked.values()))
Counter({1: 999845, 2: 69, 3: 4, 4: 1})

The marked composites take 32 MB and the marking primes take 28 MB:

>>> sum(map(getsizeof, marked))
31983680
>>> from itertools import chain
>>> sum(map(getsizeof, chain.from_iterable(marked.values())))
27999972

Together that's indeed 232 MB as reported by tracemalloc in the table above:

>>> 84 + 88 + 32 + 28
232

Since most lists of marking primes just contain a single prime anyway, let's store only the first prime at the composite. If more primes want to mark that composite, then instead move them ahead to their next multiple:

no_lists:

def primes():
    marked = {}
    for x in count(2):
        if x in marked:
            p = marked.pop(x)
            while (x := x + p) in marked:
                pass
            marked[x] = p
        else:
            yield x
            marked[x * x] = x

As the above table showed, this made it faster and saved a lot of memory.

In my other answer, I showed that prematurely setting marked[x * x] as soon as you discover a new prime x is wasteful. Let's rewrite the better version from there as an infinite generator as well. Here m is the next marking prime, and it will start marking as soon as we reach m2=m².

lazier:

def primes():
    yield 2
    primes = [2]
    markers = iter(primes)
    marked = {}
    m = next(markers)
    m2 = m * m
    for x in count(3):
        if x == m2:
            while (x := x + m) in marked:
                pass
            marked[x] = m
            m = next(markers)
            m2 = m * m
        elif x in marked:
            p = marked.pop(x)
            while (x := x + p) in marked:
                pass
            marked[x] = p
        else:
            yield x
            primes.append(x)

This brought it down to under 37 MB, as marked got much smaller. However, it creates a primes list again, in order to update m to the next next marking prime.

Can we get the best of both worlds? Not prematurely enter marked[x * x] and not build a primes list? Yes we can!

What do we use the primes list for? Just as a more slowly used stream of primes. But that's what we're writing! A generator of primes. So we can use it recursively. Our main generator iterator of primes will use a slower generator iterator of primes. But where does that slower one get its marker primes from? From yet another even more slowly progressing one, if needed. And so on. A potentially infinite stack of infinite iterators.

recursive:

(also see the addendum at the end of the answer for better versions of this)

def primes():
    yield 2
    marked = {}
    markers = primes()
    m = next(markers)
    m2 = m * m
    for x in count(3):
        if x == m2:
            while (x := x + m) in marked:
                pass
            marked[x] = m
            m = next(markers)
            m2 = m * m
        elif x in marked:
            p = marked.pop(x)
            while (x := x + p) in marked:
                pass
            marked[x] = p
        else:
            yield x

Let's see the levels by storing each iterator's marked in a global levels list by doing levels.append(marked) right after marked = {}. Then:

>>> levels = []
>>> millionth = next(islice(primes(), 999999, None))
>>> for marked in levels:
        print(len(marked), max(marked.values()))

546 3931
18 61
4 7
2 3
1 2

We see that there are five levels of iterators. The top-level iterator, the one we're using directly, has 546 marking primes, the largest being 3931. That makes sense, as the next larger prime is 3943, which starts marking at 3943², which is larger than the millionth prime:

>>> it = primes()
>>> 3931 in it and next(it)
3943
>>> 3943**2
15547249
>>> millionth
15485863

The next-lower level so far only needed to yield the 18 primes up to 61. Which is correct, as the next-larger prime 67 will only start marking at 67²=4489, well above the prime 3943 needed next by the top-level iterator. And so on, until the fifth and lowest level only needed its initial yield 2 so far.

With only a few hundred primes in all levels combined, it's no wonder this solution takes faaar less memory than our first "infinite" one, about 3000 times less, only around 100 KB. And its advantage will increase for larger numbers of primes, as the earlier solutions need O(n) memory while the recursive ones take something like O(sqrt(n)) memory (not sure how the density of primes affects it).

Last but not least, here's a little variation where I handle the prime 2 separately and then move everything twice as fast, which makes the overall solution about twice as fast (again see the table above):

twostep:

def primes():
    yield 2
    yield 3
    marked = {}
    markers = primes()
    next(markers)
    m = next(markers)
    m2 = m * m
    for x in count(5, 2):
        if x == m2:
            marked[m2 + 2*m] = 2*m
            m = next(markers)
            m2 = m * m
        elif x in marked:
            step = marked.pop(x)
            while (x := x + step) in marked:
                pass
            marked[x] = step
        else:
            yield x

Code for benchmark 1:

from timeit import repeat
from itertools import count, islice
from collections import defaultdict

def original(n):
    sieve = count(2)
    primes = []
    nprime = 0
    for _ in range(n):
        nprime = next(sieve)
        primes.append(nprime)
        sieve = filter(eval(f"lambda x: x % {nprime} != 0"), sieve)
    return primes

def improved1(n):
    primes = []
    sieve = (x for x in count(2) if all(x % p for p in primes))
    for _ in range(n):
        primes.append(next(sieve))
    return primes

def improved2(n):
    primes = []
    sieve = (x for x in count(2) if all(x % p for p in primes))
    primes += islice(sieve, n)
    return primes

def genuine(n):
    marked = defaultdict(list)
    primes = []
    for x in count(2):
        if x in marked:
            for p in marked.pop(x):
                marked[x + p].append(p)
        else:
            primes.append(x)
            if len(primes) == n:
                return primes
            marked[x * x].append(x)

funcs = original, improved1, improved2, genuine

n = 10000

expect = funcs[0](n)
for func in funcs[1:]:
    print(func(n) == expect, func.__name__)

for _ in range(3):
    for func in funcs:
        t = min(repeat(lambda: func(n), number=1))
        print('%.3f s ' % t, func.__name__)
    print()

Addendum - better recursive generators

I came up with a nicer way later, don't want to rewrite the answer using this now but want to share. This goes through the primes of the underlying generator in an outer loop and from one prime's square to the next prime's square in an inner loop. It's less code and makes the structure of the procedure clearer.

def primes():
    yield 2
    marks = {}
    i = 2
    for p in primes():
        for i in range(i+1, p*p):
            if i not in marks:
                yield i
            else:
                m = marks.pop(i)
                j = i + m
                while j in marks:
                    j += m
                marks[j] = m
        marks[p*p] = p

And a version with the paper's way of potentially storing several primes at the same common multiple:

def primes():
    yield 2
    marks = defaultdict(set)
    i = 2
    for p in primes():
        for i in range(i+1, p*p):
            if i not in marks:
                yield i
            else:
                for m in marks.pop(i):
                    marks[i+m].add(m)
        marks[p*p] = {p}

Answer 2

Back to the roots - be lazier!

(Taking this in another direction than my first answer, and it's different/long enough that I don't want to mix them.)

We can vastly improve your approach by not eagerly adding filters and instead only adding filters up to the square root of the current candidate number. And I found an in my opinion neat way to do that. Benchmark results:

                        n:   10,000   100,000  1,000,000  2,000,000 
      original_rooted_eval  0.075 s   2.036 s   59.669 s  162.764 s  
     original_rooted_asmod  0.081 s   2.128 s   56.001 s  151.860 s  
  original_rooted_default1  0.082 s   2.095 s   56.286 s  151.667 s  
  original_rooted_default2  0.075 s   1.895 s   50.088 s  134.353 s  
      original_rooted_rmod  0.041 s   1.050 s   28.101 s   75.745 s  
           improved_rooted  0.033 s   0.806 s   21.825 s   55.485 s  
            genuine_rooted  0.060 s   0.799 s   10.073 s   21.809 s  
                   genuine  0.073 s   0.952 s   12.662 s   27.009 s  

At n=10,000, this makes your original and improved versions about as fast or even up to twice as fast as my genuine sieve of Eratosthenes version. That's about a hundred times faster than they were before. But at n=1,000,000, the genuine one has already turned the table.

As the paper also points out, part of the slowness (much of it) comes from doing trial division with all primes lower than the current candidate number. It's similar to this single-number primality check:

def is_prime(x):
    return all(x % d for d in range(2, x))

We should instead only try the primes up the the square root of the candidate. If candidate \$x\$ had a divisor \$d>\sqrt{x}\$, then e=\$\frac{x}{d}<\sqrt{x}\$ would be another divisor, and we would've already found it below \$\sqrt{x}\$.

def is_prime(x):
    return x >= 2 and all(x % d for d in range(2, isqrt(x)+1))

Now how do we do that in our solutions? In your original, it means only creating the filters up to the square root of the current number. In my improved, it means not going through all primes but only to the square root. So: start filtering for 2 when you're at 4, start filtering for 3 when you're at 9, start filtering for 5 when you're at 25, etc. For this, keep track of the next prime square (i.e., 4, 9, 25, etc) and the index i of that prime in primes. Whenever you reach the next prime square, don't count it as prime but instead add the next filter.

Rooting yours:

def original_rooted_eval(n):
    sieve = count(2)
    primes = []
    next_prime_square, i = 4, 0
    for _ in range(n):
        while (nprime := next(sieve)) == next_prime_square:
            sieve = filter(eval(f"lambda x: x % {primes[i]} != 0"), sieve)
            i += 1
            next_prime_square = primes[i] ** 2
        primes.append(nprime)
    return primes

Or with your asmod:

            sieve = filter(asmod(primes[i]), sieve)

Or with a lambda with default argument:

            sieve = filter(lambda x, p=primes[i]: x % p != 0, sieve)

And btw the != 0 isn't necessary and only slows it down. We can just use the truth of the remainder:

            sieve = filter(lambda x, p=primes[i]: x % p, sieve)

And now it's obvious that we can simply use the prime's right-side modulo method:

            sieve = filter(primes[i].__rmod__, sieve)

That's all the different versions of your original that you see in the benchmark results above.

Rooting my improvement:

I could've kept for p in primes and added some if p*p > x: break, but that would be an extra cost at every innermost iteration. So instead I build an extra list of already filtering primes:

def improved_rooted(n):
    primes = []
    filters = []
    next_prime_square, i = 4, 0
    for x in count(2):
        if x == next_prime_square:
            filters.append(primes[i])
            i += 1
            next_prime_square = primes[i] ** 2
        else:
            for p in filters:
                if not x % p:
                    break
            else:
                primes.append(x)
                if len(primes) == n:
                    return primes

Rooting genuine:

My genuine already added a new prime \$x\$ into marked[x * x] so that it won't start filtering until \$x^2\$ is reached. So in that sense, genuine is already rooted. But if we never go as far as \$x^2\$, then even that is wasted. So delay putting it into marked at all until its square is reached. And then put it into its next multiple, x + primes[i] (where that x is now the square of the original, which is primes[i], so that's really primes[i]**2 + primes[i]).

def genuine_rooted(n):
    primes = []
    marked = defaultdict(list)
    next_prime_square, i = 4, 0    
    for x in count(2):
        if x == next_prime_square:
            marked[x + primes[i]].append(primes[i])
            i += 1
            next_prime_square = primes[i] ** 2
        elif x in marked:
            for p in marked.pop(x):
                marked[x + p].append(p)
        else:
            primes.append(x)
            if len(primes) == n:
                return primes

Benchmark code including full solutions:

from timeit import timeit
from itertools import count
from collections import defaultdict
from math import isqrt

def original_rooted_eval(n):
    sieve = count(2)
    primes = []
    next_prime_square, i = 4, 0
    for _ in range(n):
        while (nprime := next(sieve)) == next_prime_square:
            sieve = filter(eval(f"lambda x: x % {primes[i]} != 0"), sieve)
            i += 1
            next_prime_square = primes[i] ** 2
        primes.append(nprime)
    return primes

def original_rooted_asmod(n):
    sieve = count(2)
    primes = []
    next_prime_square, i = 4, 0
    def asmod(a):
        def inner(x):
            return x % a != 0
        return inner
    for _ in range(n):
        while (nprime := next(sieve)) == next_prime_square:
            sieve = filter(asmod(primes[i]), sieve)
            i += 1
            next_prime_square = primes[i] ** 2
        primes.append(nprime)
    return primes

def original_rooted_default1(n):
    sieve = count(2)
    primes = []
    next_prime_square, i = 4, 0
    for _ in range(n):
        while (nprime := next(sieve)) == next_prime_square:
            sieve = filter(lambda x, p=primes[i]: x % p != 0, sieve)
            i += 1
            next_prime_square = primes[i] ** 2
        primes.append(nprime)
    return primes

def original_rooted_default2(n):
    sieve = count(2)
    primes = []
    next_prime_square, i = 4, 0
    for _ in range(n):
        while (nprime := next(sieve)) == next_prime_square:
            sieve = filter(lambda x, p=primes[i]: x % p, sieve)
            i += 1
            next_prime_square = primes[i] ** 2
        primes.append(nprime)
    return primes

def original_rooted_rmod(n):
    sieve = count(2)
    primes = []
    next_prime_square, i = 4, 0
    for _ in range(n):
        while (nprime := next(sieve)) == next_prime_square:
            sieve = filter(primes[i].__rmod__, sieve)
            i += 1
            next_prime_square = primes[i] ** 2
        primes.append(nprime)
    return primes

def improved_rooted(n):
    primes = []
    filters = []
    next_prime_square, i = 4, 0
    for x in count(2):
        if x == next_prime_square:
            filters.append(primes[i])
            i += 1
            next_prime_square = primes[i] ** 2
        else:
            for p in filters:
                if not x % p:
                    break
            else:
                primes.append(x)
                if len(primes) == n:
                    return primes

def genuine_rooted(n):
    primes = []
    marked = defaultdict(list)
    next_prime_square, i = 4, 0    
    for x in count(2):
        if x == next_prime_square:
            marked[x + primes[i]].append(primes[i])
            i += 1
            next_prime_square = primes[i] ** 2
        elif x in marked:
            for p in marked.pop(x):
                marked[x + p].append(p)
        else:
            primes.append(x)
            if len(primes) == n:
                return primes

def genuine(n):
    marked = defaultdict(list)
    primes = []
    for x in count(2):
        if x in marked:
            for p in marked.pop(x):
                marked[x + p].append(p)
        else:
            primes.append(x)
            if len(primes) == n:
                return primes
            marked[x * x].append(x)

funcs = [
    original_rooted_eval,
    original_rooted_asmod,
    original_rooted_default1,
    original_rooted_default2,
    original_rooted_rmod,
    improved_rooted,
    genuine_rooted,
    genuine,
    ]

# Correctness check
n = 10**4
expect = funcs[0](n)
for func in funcs[1:]:
    print(func(n) == expect, func.__name__)
print()

width = max(len(func.__name__) for func in funcs) + 2

sizes = 10**4, 10**5, 10**6, 2 * 10**6
rounds = 5

# tss[s][f] := list of times for size #s and function #f.
tsss = []
for n in sizes:
    tss = [[] for _ in funcs]
    tsss.append(tss)
    for r in range(rounds):
        print(f'{n=:,}  round={r+1}/{rounds}')
        for func, ts in zip(funcs, tss):
            # t = n / 1e5  # fake for faster development
            t = timeit(lambda: func(n), number=1)
            ts.append(t)
            # print(func.__name__, t)  # Show results as they happen
        print('n:'.rjust(width), *(f'{n:8,} ' for n in sizes))
        for f, func in enumerate(funcs):
            print(func.__name__.rjust(width), end=' ')
            for tss in tsss:
                print(' %.3f s  ' % min(tss[f]), end='')
            print()
        print()

Answer 3

Bonus points for coming up with something that is both "interesting" and "technically correct"; but that's about where it stops. It's certainly not fast.

Think about what's happening here: for every single call of eval or asmod, you're chaining another layer of filter - with a new function - on your generator. That simply won't scale. You're better off with the simple approach, having a single data structure that stores your filter. Experiment with a list or array of booleans, or set of integers, to see which one has the best performance.

Answer 4

"Purist" in me says that any solution which is not starting with building finite table of numbers is not sieve of Eratosthenes, as well as any solutions which don't keep all information (numbers, cross-offs) in one table (datastructure).

Why? Back in Eratosthenes days he had to create table, after he finished creating table he started crossing off numbers (not deleting them from table). End result was the same table except non-primes (while present) were 'crossed off'.

Keeping this interpretation in mind following is more akin representation of operations Eratosthenes made:

• create finite 'table'. For that abstraction one can use np.array. It will be one row one column table where array values represent 'cross off' and indices represent numbers. In initial state all numbers are not crossed off i.e. True:

sieve = np.ones(max_value + 1, dtype=np.bool)

• 'cross off' numbers until first prime (2) by setting values to 'False'

sieve[:2] = False

• go through all numbers starting from first prime (2) up to integer part of square root + 1, if number is not 'crossed off' (i.e value of corresponding index is True), 'cross off' numbers starting of square of number with step of number.

for i in range(2, isqrt(max_value)+1): 
    if sieve[i]:
        sieve[i*i::i] = False

We have 'table' where indices which are 'crossed off' have value False and all 'non-crossed off' indices values are True.

For me this way is more mimicking paper-and-pencil (or clay and stylus) counterpart of original sieve of Eratosthenes (your mileage may vary). As simple function:

import numpy as np
from math import isqrt


def primes(max_value):
    sieve = np.ones(max_value + 1, dtype=np.bool)
    sieve[:2] = False
    for i in range(2, isqrt(max_value)+1): 
        if sieve[i]:
            sieve[i*i::i] = False
    return np.nonzero(sieve)[0]            # return only primes

And now something completely different: on holy page of Life, Universe and Everything else (python.org) there is SimplePrograms: "20 lines: Prime numbers sieve w/fancy generators":

import itertools

def iter_primes():
     # an iterator of all numbers between 2 and +infinity
     numbers = itertools.count(2)

     # generate primes forever
     while True:
         # get the first number from the iterator (always a prime)
         prime = next(numbers)
         yield prime

         # this code iteratively builds up a chain of
         # filters...slightly tricky, but ponder it a bit
         numbers = filter(prime.__rmod__, numbers)

for p in iter_primes():
    if p > 1000:
        break
    print (p)