An Efficient Sieve of Eratosthenes in Python

Question

This very short and simple code in #Python tries to simulate the "Sieve of Eratosthenes" for the first N natural numbers with the constraints of (0) script shortness; (1) minimization of the 'if statements' and 'for/while loops'; (2) efficiency in terms of CPU time.

import numpy as np
N = 10**5
a = np.array(range(3,N,2))
for j in range(0, int(round(np.sqrt(N),0))):
    a[(a!=a[j]) & (a%a[j] == 0)] = 0
    a = a[a!=0]
a = [2]+list(a)

On an Intel Core I5, it returns the prime numbers among the first:

• N = 100,000 in 0.03 seconds;
• N = 1,000,000 in 0.63 seconds;
• N = 10,000,000 in 22.2 seconds.

Would someone like to share more efficient codes in term of CPU time within the aforementioned constraints?

Answer 1

An actual NumPy sieve of Eratosthenes looks like this:

def sieve(n):
    flags = numpy.ones(n, dtype=bool)
    flags[0] = flags[1] = False
    for i in range(2, n):
        # We could use a lower upper bound for this loop, but I don't want to bother with
        # getting the rounding right on the sqrt handling.
        if flags[i]:
            flags[i*i::i] = False
    return numpy.flatnonzero(flags)

It maintains an array of "possibly prime" flags and directly unsets the flags corresponding to multiples of primes, without needing to test divisibility, especially for numbers that aren't divisible by the prime currently being handled.

What you're doing is trial division, where you just go through and test whether numbers are divisible by candidate divisors. Even a good implementation of trial division needs to do more operations, and more expensive operations, than a sieve. Your implementation does even more work than that, because it considers non-prime candidate divisors, and because it keeps performing divisibility tests for numbers it should already know are prime.