Why is my Sieve of Eratosthenes using generators so slow

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.

• This is not the sieve of Eratosthenes Oct 26, 2020 at 17:22
• @trentcl: One of my favorite papers, which should be much more widely known. Oct 27, 2020 at 5:42
• You could still use the lambda without relying on eval; you just need to prebind nprime as a defaulted argument: sieve = filter(lambda x, nprime=nprime: x % nprime != 0, sieve). That said, if you need a Python defined function (def or lambda) to use filter, you're better off just using an equivalent listcomp or genexpr (sadly, harder here due to late binding), or finding a way to make it actually use a built-in implemented in C. In this case, you could just do sieve = filter(nprime.__rmod__, sieve) with the same effect (given nprime and all sieve elements are int). Oct 27, 2020 at 18:42
• Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see *what you may and may not do after receiving answers. Oct 28, 2020 at 16:43
• @EricDuminil I guess the Python version would use itertools to do wheel2357 = cycle([2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10]) and primes = chain([2, 3, 5, 7], sieve(accumulate(wheel2357, initial=11))).. Nov 14, 2020 at 0:01

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=m2.

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 39432, 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 672=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[p*p] = {p}

• To be clear, the first improvement is converting the algorithm from recursion to iteration, avoiding the overhead of calling all the lambdas, the actual test operations are the same. I feel like the phrasing isn't super clear since "lets 2/3 of all numbers through" can suggest that 6 would be checked by filter_2 AND filter_3 rather than being ejected at the filter_2 test. I had to re-read a couple of times to twig that you describe each filter in isolation then show the result of stacking them. Really, it's more of a grammar complaint than anything. Oct 28, 2020 at 14:03
• That's what I mean by changing it from recursive to iterative. The way I double checked the percents was the opposite way... since 1/3rd of multiples of 2 are multiples of 3 (obviously), the proportion of multiples in the removed and remaining is the same, thus each filter is independent of its predecessors if the input set is uniform. But like I said, it's purely the pedantic difference between implicit "tests all numbers" and explicit "tests remaining numbers". What's really important is how elegant the marking-next-square approach to the genuine really is. Oct 28, 2020 at 16:23
• As a bonus, you can easily convert it to an infinite prime generator : tio.run/##bY/NasMwEITveoqBYiI1JiTNLZB36L2UYqx1u1SWFlkBpS/… Oct 28, 2020 at 18:30
• @EricDuminil I didn't read that part of the paper, partly because I don't understand Haskell so at some point didn't understand anymore what they're saying. I did try a mini-wheel with just 2 and 3 and I considered a larger wheel, but then felt like the solutions grew complicated enough already. That said, I found another improvement that made it a bit faster and simpler, so maybe I'll try to incorporate a wheel with that. Nov 1, 2020 at 11:32
• @EricDuminil Didn't manage to incorporate larger wheels, but added the improved version(s) at the end now. Mar 9, 2021 at 22:00

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()


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.

• How would you write it with list/array/set? Oct 26, 2020 at 17:13
• @HeapOverflow That's such a common question that there's an entire tag dedicated to it: codereview.stackexchange.com/questions/tagged/… Oct 26, 2020 at 17:15
• But they're not really doing the sieve of Eratosthenes. I was wondering about your version of their fake one. Oct 26, 2020 at 17:30
• @Reinderien I was mostly thinking to make a generator that can go indefinitely, which would not work with a list, array, etc. Oct 26, 2020 at 17:38
• @jayjayjay It's possible to have a segmented, progressive sieve that can go indefinitely but is still data-structure-based. Oct 26, 2020 at 17:55

"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)