Sieve of Eratosthenes#
Programming Environment#
import math
from numba import jit, njit, prange
import numpy as np
from typing import Iterator
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 2
1 import math
----> 2 from numba import jit, njit, prange
3 import numpy as np
4 from typing import Iterator
ModuleNotFoundError: No module named 'numba'
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)
25
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.4 µs ± 57.9 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
204 µs ± 398 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
929 ms ± 6.23 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
n = int(1e8)
%time sieve = [True] * n
%time sieve = [True for _ in range(n)]
CPU times: user 19.7 ms, sys: 54.4 ms, total: 74 ms
Wall time: 73.7 ms
CPU times: user 1.34 s, sys: 73.2 ms, total: 1.41 s
Wall time: 1.41 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}\]
# 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]])
3 84 [ 0 1 2 3 4 5 6 7 8 10 11 12 13 14 16 17 18 19 20 22 23 24 25 26
28 29 30 31 32 34 35 36 37 38 40 41 42 43 44 46 47 48 49 50 52 53 54 55
56 58 59 60 61 62 64 65 66 67 68 70 71 72 73 74 76 77 78 79 80 82 83 84
85 86 88 89 90 91 92 94 95 96 97 98]
5 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
29 30 31 32 34 36 37 38 40 41 42 43 44 46 47 48 49 50 52 53 54 56 58 59
60 61 62 64 66 67 68 70 71 72 73 74 76 77 78 79 80 82 83 84 86 88 89 90
91 92 94 96 97 98]
7 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
29 30 31 32 34 36 37 38 40 41 42 43 44 46 47 48 50 52 53 54 56 58 59 60
61 62 64 66 67 68 70 71 72 73 74 76 78 79 80 82 83 84 86 88 89 90 92 94
96 97 98]
9 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
29 30 31 32 34 36 37 38 40 41 42 43 44 46 47 48 50 52 53 54 56 58 59 60
61 62 64 66 67 68 70 71 72 73 74 76 78 79 80 82 83 84 86 88 89 90 92 94
96 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))
197 ms ± 2.72 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
193 ms ± 2.2 ms per loop (mean ± std. dev. of 7 runs, 10 loops 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))
11.5 ms ± 44.3 µ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))
---------------------------------------------------------------------------
UnsupportedRewriteError Traceback (most recent call last)
/Users/df/Desktop/._characteristica_universalis/my_github/12___courses_public/logic-and-numbers/numbers/src/sieve_of_eratosthenes.ipynb Cell 29 line 1
<a href='vscode-notebook-cell:/Users/df/Desktop/._characteristica_universalis/my_github/12___courses_public/logic-and-numbers/numbers/src/sieve_of_eratosthenes.ipynb#X41sZmlsZQ%3D%3D?line=7'>8</a> primes[i*i::i] = False
<a href='vscode-notebook-cell:/Users/df/Desktop/._characteristica_universalis/my_github/12___courses_public/logic-and-numbers/numbers/src/sieve_of_eratosthenes.ipynb#X41sZmlsZQ%3D%3D?line=8'>9</a> return np.flatnonzero(a = primes)
---> <a href='vscode-notebook-cell:/Users/df/Desktop/._characteristica_universalis/my_github/12___courses_public/logic-and-numbers/numbers/src/sieve_of_eratosthenes.ipynb#X41sZmlsZQ%3D%3D?line=10'>11</a> get_ipython().run_line_magic('timeit', 'p5(int(1e7))')
File ~/anaconda3/envs/ml/lib/python3.11/site-packages/IPython/core/interactiveshell.py:2432, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2430 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2431 with self.builtin_trap:
-> 2432 result = fn(*args, **kwargs)
2434 # The code below prevents the output from being displayed
2435 # when using magics with decorator @output_can_be_silenced
2436 # when the last Python token in the expression is a ';'.
2437 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File ~/anaconda3/envs/ml/lib/python3.11/site-packages/IPython/core/magics/execution.py:1185, in ExecutionMagics.timeit(self, line, cell, local_ns)
1183 for index in range(0, 10):
1184 number = 10 ** index
-> 1185 time_number = timer.timeit(number)
1186 if time_number >= 0.2:
1187 break
File ~/anaconda3/envs/ml/lib/python3.11/site-packages/IPython/core/magics/execution.py:173, in Timer.timeit(self, number)
171 gc.disable()
172 try:
--> 173 timing = self.inner(it, self.timer)
174 finally:
175 if gcold:
File <magic-timeit>:1, in inner(_it, _timer)
File ~/anaconda3/envs/ml/lib/python3.11/site-packages/numba/core/dispatcher.py:471, in _DispatcherBase._compile_for_args(self, *args, **kws)
468 error_rewrite(e, 'typing')
469 except errors.UnsupportedError as e:
470 # Something unsupported is present in the user code, add help info
--> 471 error_rewrite(e, 'unsupported_error')
472 except (errors.NotDefinedError, errors.RedefinedError,
473 errors.VerificationError) as e:
474 # These errors are probably from an issue with either the code
475 # supplied being syntactically or otherwise invalid
476 error_rewrite(e, 'interpreter')
File ~/anaconda3/envs/ml/lib/python3.11/site-packages/numba/core/dispatcher.py:409, in _DispatcherBase._compile_for_args.<locals>.error_rewrite(e, issue_type)
407 raise e
408 else:
--> 409 raise e.with_traceback(None)
UnsupportedRewriteError: Failed in nopython mode pipeline (step: convert to parfors)
Only constant step size of 1 is supported for prange
File "../../../../../../../../../var/folders/89/5r24znsj4jbfr7zrccy5yn0c0000gn/T/ipykernel_68182/4269974924.py", line 6:
<source missing, REPL/exec in use?>
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))