Sieve of Eratosthenes - Finding Primes Python

Viewed 120977

Just to clarify, this is not a homework problem :)

I wanted to find primes for a math application I am building & came across Sieve of Eratosthenes approach.

I have written an implementation of it in Python. But it's terribly slow. For say, if I want to find all primes less than 2 million. It takes > 20 mins. (I stopped it at this point). How can I speed this up?

def primes_sieve(limit):
    limitn = limit+1
    primes = range(2, limitn)

    for i in primes:
        factors = range(i, limitn, i)
        for f in factors[1:]:
            if f in primes:
                primes.remove(f)
    return primes

print primes_sieve(2000)

UPDATE: I ended up doing profiling on this code & found that quite a lot of time was spent on removing an element from the list. Quite understandable considering it has to traverse the entire list (worst-case) to find the element & then remove it and then readjust the list (maybe some copy goes on?). Anyway, I chucked out list for dictionary. My new implementation -

def primes_sieve1(limit):
    limitn = limit+1
    primes = dict()
    for i in range(2, limitn): primes[i] = True

    for i in primes:
        factors = range(i,limitn, i)
        for f in factors[1:]:
            primes[f] = False
    return [i for i in primes if primes[i]==True]

print primes_sieve1(2000000)
22 Answers

Using a bit of numpy, I could find all primes below 100 million in a little over 2 seconds.

There are two key features one should note

  • Cut out multiples of i only for i up to root of n
  • Setting multiples of i to False using x[2*i::i] = False is much faster than an explicit python for loop.

These two significantly speed up your code. For limits below one million, there is no perceptible running time.

import numpy as np

def primes(n):
    x = np.ones((n+1,), dtype=np.bool)
    x[0] = False
    x[1] = False
    for i in range(2, int(n**0.5)+1):
        if x[i]:
            x[2*i::i] = False

    primes = np.where(x == True)[0]
    return primes

print(len(primes(100_000_000)))
import math
def sieve(n):
    primes = [True]*n
    primes[0] = False
    primes[1] = False
    for i in range(2,int(math.sqrt(n))+1):
            j = i*i
            while j < n:
                    primes[j] = False
                    j = j+i
    return [x for x in range(n) if primes[x] == True]

The fastest implementation I could come up with:

isprime = [True]*N
isprime[0] = isprime[1] = False
for i in range(4, N, 2):
    isprime[i] = False
for i in range(3, N, 2):
    if isprime[i]:
        for j in range(i*i, N, 2*i):
            isprime[j] = False

i think this is shortest code for finding primes with eratosthenes method

def prime(r):
    n = range(2,r)
    while len(n)>0:
        yield n[0]
        n = [x for x in n if x not in range(n[0],r,n[0])]


print(list(prime(r)))

I just came up with this. It may not be the fastest, but I'm not using anything other than straight additions and comparisons. Of course, what stops you here is the recursion limit.

def nondivsby2():
    j = 1
    while True:
        j += 2
        yield j

def nondivsbyk(k, nondivs):
    j = 0
    for i in nondivs:
        while j < i:
            j += k
        if j > i:
            yield i

def primes():
    nd = nondivsby2()
    while True:
        p = next(nd)
        nd = nondivsbyk(p, nd)
        yield p

def main():
    for p in primes():
        print(p)

I made a one liner version of the Sieve of Eratosthenes

sieve = lambda j: [print(x) for x in filter(lambda n: 0 not in map(lambda i: n % i, range(2, n)) and (n!=1)&(n!=0), range(j + 1))]

In terms of performance, I am pretty sure this isn't the fastest thing by any means, and in terms of readability / following PEP8, this is pretty terrible, but it's more the novelty of the length than anything.

EDIT: Note that this simply prints the sieve & does not return (if you attempt to print it you will get a list of Nones, if you want to return, change the print(x) in the list comprehension to just "x".

not sure if my code is efficeient, anyone care to comment?

from math import isqrt

def isPrime(n):
    if n >= 2: # cheating the 2, is 2 even prime?
        for i in range(3, int(n / 2 + 1),2): # dont waste time with even numbers
            if n % i == 0:
                return False
    return True

def primesTo(n): 
    x = [2] if n >= 2 else [] # cheat the only even prime
    if n >= 2:
        for i in range(3, n + 1,2): # dont waste time with even numbers
            if isPrime(i):
                x.append(i)  
    return x

def primes2(n): # trying to do this using set methods and the "Sieve of Eratosthenes"
    base = {2} # again cheating the 2
    base.update(set(range(3, n + 1, 2))) # build the base of odd numbers
    for i in range(3, isqrt(n) + 1, 2): # apply the sieve
        base.difference_update(set(range(2 * i, n + 1 , i)))
    return list(base)

print(primesTo(10000)) # 2 different methods for comparison
print(primes2(10000))

Probably the quickest way to have primary numbers is the following:

import sympy
list(sympy.primerange(lower, upper+1))

In case you don't need to store them, just use the code above without conversion to the list. sympy.primerange is a generator, so it does not consume memory.

Using recursion and walrus operator:

def prime_factors(n):
    for i in range(2, int(n ** 0.5) + 1):
        if (q_r := divmod(n, i))[1] == 0:
            return [i] + factor_list(q_r[0])
    return [n]

Basic sieve

with numpy is amazing fast. May be the fastest implementation

# record: sieve 1_000_000_000 in 6.9s (core i7 - 2.6Ghz)
def sieve_22max_naive(bound):
    sieve = np.ones(bound, dtype=bool)  # default all prime
    sieve[:2] = False  # 0, 1 is not prime

    sqrt_bound = math.ceil(math.sqrt(bound))

    for i in range(2, sqrt_bound):
        if sieve[i]:
            inc = i if i == 2 else 2 * i
            sieve[i * i:bound:inc] = False

    return np.arange(bound)[sieve]


if __name__ == '__main__':
    start = time.time()
    prime_list = sieve_22max_naive(1_000_000_000)
    print(f'Count: {len(prime_list):,}\n'
          f'Greatest: {prime_list[-1]:,}\n'
          f'Elapsed: %.3f' % (time.time() - start))

Segment sieve (use less memory)

# find prime in range [from..N), base on primes in range [2..from)
def sieve_era_part(primes, nfrom, n):
    sieve_part = np.ones(n - nfrom, dtype=bool)  # default all prime

    limit = math.ceil(math.sqrt(n))
    # [2,3,5,7,11...p] can find primes < (p+2)^2
    if primes[-1] < limit - 2:
        print(f'Not enough base primes to find up to {n:,}')
        return

    for p in primes:
        if p >= limit: break

        mul = p * p
        inc = p * (2 if p > 2 else 1)
        if mul < nfrom:
            mul = math.ceil(nfrom / p) * p
            (mul := mul + p) if p > 2 and (mul & 1) == 0 else ...  # odd, not even

        sieve_part[mul - nfrom::inc] = False

    return np.arange(nfrom, n)[sieve_part]
    # return np.where(sieve_part)[0] + nfrom
    # return [i + nfrom for i, is_p in enumerate(sieve_part) if is_p]
    # return [i for i in range(max(nfrom, 2), n) if sieve_part[i - nfrom]]


# find nth prime number, use less memory,
# extend bound to SEG_SIZE each loop
# record: 50_847_534 nth prime in 6.78s, core i7 - 9850H 2.6GHhz
def nth_prime(n):
    # find prime up to bound
    bound = 500_000
    primes = sieve_22max_naive(bound)

    SEG_SIZE = int(50e6)

    while len(primes) < n:
        # sieve for next segment
        new_primes = sieve_era_part(primes, bound, bound + SEG_SIZE)
        # extend primes
        bound += SEG_SIZE
        primes = np.append(primes, new_primes)

    return primes[n - 1]


if __name__ == '__main__':
    start = time.time()
    prime = nth_prime(50_847_534)
    print(f'{prime:,} Time %.6f' % (time.time() - start))

here is my solution, the same as Wikipedia

import math


def sieve_of_eratosthenes(n):
    a = [i for i in range(2, n+1)]
    clone_a = a[:]
    b = [i for i in range(2, int(math.sqrt(n))+1)]
    for i in b:
        if i in a:
            c = [pow(i, 2)+(j*i) for j in range(0, n+1)]
            for j in c:
                if j in clone_a:
                    clone_a.remove(j)
    return clone_a


if __name__ == '__main__':
    print(sieve_of_eratosthenes(23))
Related