Sieve of Eratosthenes

Sieve of Eratosthenes#

https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n/2068548

Programming Environment#

import math
#from   numba  import jit, njit, prange
import numpy  as np
from   typing import Iterator

def generate_integers (start : int = 0) -> Iterator[int]:
  n = start
  while True:
    yield n
    n += 1

num = generate_integers()
print(list(next(num) for _ in range(25)))

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]

next(num)

Naive Python#

def p0 (p    : int,
        flag : bool = False) -> list[int]:
  num = generate_integers(2)
  arr = list(next(num) for _ in range(2, p + 1))

  primes = []
  prime  = arr.pop(0)  # remove the first number in the array of numbers
  primes.append(prime) # add it to our list of primes
  
  while prime <= p ** 0.5:
    if flag:
      print(prime, arr)
    arr   = list(i for i in arr if i % primes[-1] != 0) # remove all prime multiples from the array:
    prime = arr.pop(0)
    primes.append(prime)

  primes += arr
  return primes

%timeit p0(int(1e2))
%timeit p0(int(1e3))
%timeit p0(int(1e6))

15.6 µs ± 62.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

210 µs ± 1.19 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

1.02 s ± 20.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

p0(int(1e6))

n = int(1e8)
%time sieve = [True] * n
%time sieve = [True for _ in range(n)]

CPU times: user 20.1 ms, sys: 50.7 ms, total: 70.8 ms
Wall time: 71.9 ms

CPU times: user 1.3 s, sys: 127 ms, total: 1.43 s
Wall time: 1.46 s

# def sieve_of_eratosthenes (n):
#   sieve      = [True] * n
#   sieve[0:1] = [False, False]
#   for prime in range(2, n + 1):
#     if sieve[prime]:
#       for i in range(2*prime, n+1, prime):
#         sieve[i] = False

# def p1 (n):
#   sieve = [True] * n
#   p     = 2
#   while (p*p <= n):
#     if (sieve[p] == True):
#       for i in range(p*p, n+1, p):
#         sieve[i] = False
#     p += 1
#   for p in range(2, n+1):
#     if sieve[p]:
#       p

# %timeit p1(int(1e2))
# %timeit p1(int(1e3))
# %timeit p1(int(1e6))
# %timeit p1(int(1e7))

\[\begin{split} \begin{align*} n = 3 && 9 + 6k = 9, 15, 21, 27, 33, 39, 45, 51, 57, 63, 69, 75, 81, 87, 93, 99 \\ n = 5 && 25 + 10k = 25, 35, 45, 55, 65, 75, 85, 95 \\ n = 7 && 49 + 14k = 49, 63, 77, 91 \\ n = 11 && 121 + 22k = \dots \\ \end{align*} \end{split}\]

\( \begin{aligned} & 3^2 + 2.3k = 3(2k + 3) = \phantom{1} 6k + \phantom{1} 9 = \phantom{1} 6(k + 1) + \phantom{1} 3 = 3(2k' + 1) && \text{odd multiples of 3} \\ & 5^2 + 2.5k = 5(2k + 5) = 10k + 25 = 10(k + 2) + \phantom{1} 5 = 5(2k' + 1) && \text{odd multiples of 5 but not 3} \\ & 7^2 + 2.7k = 7(2k + 7) = 14k + 49 = 14(k + 3) + \phantom{1} 7 = 7(2k' + 1) && \text{odd multiples of 7 but not 3, 5} \\ \end {aligned} \)

# FULL SIEVE
n     = 100
sieve = [True] * n
for i in range(3, int(n ** 0.5) + 1, 2):
  if sieve[i]:
    sieve[i*i::2*i] = [False] * ((n - i * i - 1) // (2 * i) + 1)
  print(i, len(np.flatnonzero(sieve)), np.flatnonzero(sieve))
print([2] + [i for i in range(3, n, 2) if sieve[i]])

[ 0  1  2  3  4  5  6  7  8 10 11 12 13 14 16 17 18 19 20 22 23 24 25 26
29 30 31 32 34 35 36 37 38 40 41 42 43 44 46 47 48 49 50 52 53 54 55
58 59 60 61 62 64 65 66 67 68 70 71 72 73 74 76 77 78 79 80 82 83 84
86 88 89 90 91 92 94 95 96 97 98]
78 [ 0  1  2  3  4  5  6  7  8 10 11 12 13 14 16 17 18 19 20 22 23 24 26 28
30 31 32 34 36 37 38 40 41 42 43 44 46 47 48 49 50 52 53 54 56 58 59
61 62 64 66 67 68 70 71 72 73 74 76 77 78 79 80 82 83 84 86 88 89 90
92 94 96 97 98]
75 [ 0  1  2  3  4  5  6  7  8 10 11 12 13 14 16 17 18 19 20 22 23 24 26 28
30 31 32 34 36 37 38 40 41 42 43 44 46 47 48 50 52 53 54 56 58 59 60
62 64 66 67 68 70 71 72 73 74 76 78 79 80 82 83 84 86 88 89 90 92 94
97 98]
75 [ 0  1  2  3  4  5  6  7  8 10 11 12 13 14 16 17 18 19 20 22 23 24 26 28
30 31 32 34 36 37 38 40 41 42 43 44 46 47 48 50 52 53 54 56 58 59 60
62 64 66 67 68 70 71 72 73 74 76 78 79 80 82 83 84 86 88 89 90 92 94
97 98]
[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]

# HALF SIEVE
n     = 100
sieve = [True] * (n//2)
for i in range(3, int(n**0.5) + 1, 2):
  if sieve[i//2]:
    sieve[i*i//2::i] = [False] * ((n - i*i - 1) // (2*i) + 1)
  print(i, len(np.flatnonzero(sieve)), np.flatnonzero(sieve))
print([2] + [2*i + 1 for i in range(1, n//2) if sieve[i]])

3 34 [ 0  1  2  3  5  6  8  9 11 12 14 15 17 18 20 21 23 24 26 27 29 30 32 33
 35 36 38 39 41 42 44 45 47 48]
5 28 [ 0  1  2  3  5  6  8  9 11 14 15 18 20 21 23 24 26 29 30 33 35 36 38 39
 41 44 45 48]
7 25 [ 0  1  2  3  5  6  8  9 11 14 15 18 20 21 23 26 29 30 33 35 36 39 41 44
 48]
9 25 [ 0  1  2  3  5  6  8  9 11 14 15 18 20 21 23 26 29 30 33 35 36 39 41 44
 48]
[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]

def p2 (n):
  """ FULL SIEVE
  Returns a list of primes p < n.
  https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n/2068548
  """
  sieve = [True] * n
  for i in range(3, int(n ** 0.5) + 1, 2):
    if sieve[i]:
      sieve[i*i::2*i] = [False] * ((n - i * i - 1) // (2 * i) + 1)
  return [2] + [i for i in range(3, n, 2) if sieve[i]]

def p3 (n):
  """ HALF SIEVE
  Returns a list of primes p < n.
  https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n/2068548
  """
  sieve = [True] * (n//2)
  for i in range(3, int(n**0.5) + 1, 2):
    if sieve[i//2]:
      sieve[i*i//2::i] = [False] * ((n - i*i - 1) // (2*i) + 1)
  return [2] + [2*i + 1 for i in range(1, n//2) if sieve[i]]

%timeit p2(int(1e7))
%timeit p3(int(1e7))

208 ms ± 5.13 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

209 ms ± 3.13 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

NumPy and Numba#

def p4 (n : int) -> np.ndarray:
  """Sieve of Eratosthenes"""
  flags = np.full(shape = n, fill_value = True) #np.ones(shape = n, dtype = bool)
  flags[0], flags[1] = False, False # 0 and 1 are not prime
  for i in range(2, int(n ** 0.5) + 1, 2):
    if flags[i]:
      flags[i*i::i] = False
  return np.flatnonzero(a = flags)

%timeit p4(int(1e7))

12.7 ms ± 89.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

@njit(parallel = True, fastmath = True)
def p5 (n : int) -> np.ndarray:
  """Sieve of Eratosthenes"""
  primes = np.full(shape = n, fill_value = True)
  primes[0], primes[1] = False, False # 0 and 1 are not prime
  for i in prange(2, int(n ** 0.5) + 1, 2):
    if primes[i]:
      primes[i*i::i] = False
  return np.flatnonzero(a = primes)

%timeit p5(int(1e7))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[17], line 1
----> 1 @njit(parallel = True, fastmath = True)
      2 def p5 (n : int) -> np.ndarray:
      3   """Sieve of Eratosthenes"""
      4   primes = np.full(shape = n, fill_value = True)

NameError: name 'njit' is not defined

def p6 (n : int) -> np.ndarray:
  """ Returns an array of primes, 3 <= p < n """
  sieve = np.full(shape = n//2, fill_value = True) #np.ones(n//2, dtype=bool)
  for i in range(3, int(n**0.5) + 1, 2):
    if sieve[i//2]:
      sieve[i*i//2::i] = False
  return 2*np.nonzero(a = sieve)[0][1::] + 1

%timeit p6(int(1e7))

2.5 s ± 67.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

def p7 (n : int) -> np.ndarray:
  """ Input n>=6, Returns an array of primes, 2 <= p < n """
  sieve = np.full(shape = n//3 + (n%6==2), fill_value = True) #np.ones(n//3 + (n%6==2), dtype=bool)
  for i in range(3, int(n**0.5)//3 + 1):
    if sieve[i]:
      k = 3*i + 1 | 1
      sieve[k*k//3            ::2*k] = False
      sieve[k*(k-2*(i&1)+4)//3::2*k] = False
  return np.r_[2,3,((3*np.nonzero(sieve)[0][1:]+1)|1)]

%timeit p7(int(1e7))

2.57 s ± 78.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit p0(int(1e7))
%timeit p1(int(1e7))
%timeit p2(int(1e7))
%timeit p3(int(1e7))
%timeit p4(int(1e7))
%timeit p5(int(1e7))
%timeit p6(int(1e7))
%timeit p7(int(1e7))

Sieve of Eratosthenes

Contents

Sieve of Eratosthenes#

Programming Environment#

Naive Python#

NumPy and Numba#