Presumo que todos vocês saibam o que são números primos e o que é o Crivo de Eratóstenes, então não vou perder tempo explicando-os.
Agora, todos os números primos, exceto 2, são números ímpares, então precisamos verificar apenas os números ímpares. Isso é muito óbvio, mas vale a pena mencionar porque essa otimização simples reduziu pela metade os candidatos que precisamos verificar, e então precisamos marcar apenas os múltiplos ímpares de números primos diferentes de 2.
Outra otimização simples é verificar apenas números até a raiz quadrada do limite, e outra otimização simples é marcar como números compostos o início no quadrado do número primo. Como os motivos para isso devem ser óbvios, não explicarei o porquê, embora eu não consiga pensar em outras otimizações relacionadas à marcação de números compostos.
Mas podemos restringir ainda mais o espaço de busca: todos os números primos, exceto 2, 3 e 5, devem ser congruentes a [1, 7, 11, 13, 17, 19, 23, 29]30%. Isso é evidente pela natureza do módulo. Existem apenas 30 possibilidades; se o módulo for par, o número é múltiplo de 2, e as outras possibilidades significam que o número é múltiplo de 3, múltiplo de 5 ou ambos. Em outras palavras, todos os números primos devem ser coprimos a 30, exceto 2, 3 e 5.
Em Python, isso é:
wheel3 = [i for i in range(1, 30) if math.gcd(i, 30) == 1]
# [1, 7, 11, 13, 17, 19, 23, 29]
Agora calculamos a diferença entre pares consecutivos de elementos, começando em 7, e 1 vem imediatamente depois de 29 devido à natureza da operação de módulo, por exemplo, 31 % 30 == 1e então a diferença entre eles é 2.
Obtemos o seguinte: [4, 2, 4, 2, 4, 6, 2, 6].
Agora, de cada 30 números, precisamos verificar apenas 8, pulando 22. Isso representa uma melhoria significativa em relação à otimização anterior, que aplicava força bruta apenas aos números ímpares. Precisávamos processar 15 números de cada 30 números se processássemos todos os números ímpares. Agora, temos 7 números a menos, e o espaço de busca é reduzido para 4/15, que é 0,2666...
Podemos otimizar isso ainda mais usando uma roda maior, usando 2 * 3 * 5 * 7 = 210 como base, todos os números primos começando em 11 devem ser coprimos com 210.
wheel4 = [i for i in range(1, 210) if math.gcd(i, 210) == 1]
wheel4 == [
1 , 11 , 13 , 17 , 19 , 23 ,
29 , 31 , 37 , 41 , 43 , 47 ,
53 , 59 , 61 , 67 , 71 , 73 ,
79 , 83 , 89 , 97 , 101, 103,
107, 109, 113, 121, 127, 131,
137, 139, 143, 149, 151, 157,
163, 167, 169, 173, 179, 181,
187, 191, 193, 197, 199, 209
]
E a lista de alterações de índice é:
[
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
]
Agora, de cada 210 números, precisamos processar apenas 48 números, em comparação com os 56 números anteriores, precisamos processar 8 números a menos, reduzindo o espaço de busca para 8/35, que é 0,22857142857142857... e menos de um quarto.
Então, espero que a versão que usa a roda baseada em 210 leve apenas 6/7 ou 85,71% do tempo que a versão baseada em 30 leva para ser executada, mas não é isso que acontece.
import math
import numpy as np
import numba as nb
@nb.njit(cache=True)
def prime_wheel_sieve(n: int) -> np.ndarray:
wheel = [4, 2, 4, 2, 4, 6, 2, 6]
primes = np.ones(n + 1, dtype=np.uint8)
primes[:2] = False
for square, step in ((4, 2), (9, 6), (25, 10)):
primes[square::step] = False
k = 7
lim = int(math.sqrt(n) + 0.5)
i = 0
while k <= lim:
if primes[k]:
primes[k**2 :: 2 * k] = False
k += wheel[i]
i = (i + 1) & 7
return np.nonzero(primes)[0]
# fmt: off
WHEEL4 = np.array([
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
], dtype=np.uint8)
# fmt: on
@nb.njit(cache=True)
def prime_wheel_sieve4(n: int) -> np.ndarray:
primes = np.ones(n + 1, dtype=np.uint8)
primes[:2] = False
for square, step in ((4, 2), (9, 6), (25, 10), (49, 14)):
primes[square::step] = False
k = 11
lim = int(math.sqrt(n) + 0.5)
i = 0
while k <= lim:
if primes[k]:
primes[k**2 :: 2 * k] = False
k += WHEEL4[i]
i = (i + 1) % 48
return np.nonzero(primes)[0]
In [549]: np.array_equal(prime_wheel_sieve(65536), prime_wheel_sieve4(65536))
Out[549]: True
In [550]: %timeit prime_wheel_sieve(65536)
161 μs ± 1.13 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [551]: %timeit prime_wheel_sieve4(65536)
163 μs ± 1.79 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [552]: %timeit prime_wheel_sieve4(131072)
330 μs ± 10.6 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [553]: %timeit prime_wheel_sieve(131072)
328 μs ± 7.4 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [554]: %timeit prime_wheel_sieve4(262144)
680 μs ± 14.3 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [555]: %timeit prime_wheel_sieve(262144)
669 μs ± 7.79 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [556]: %timeit prime_wheel_sieve(524288)
1.44 ms ± 16.2 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [557]: %timeit prime_wheel_sieve4(524288)
1.48 ms ± 13.4 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [558]: %timeit prime_wheel_sieve4(1048576)
3.25 ms ± 81.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [559]: %timeit prime_wheel_sieve(1048576)
3.23 ms ± 115 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [560]: %timeit prime_wheel_sieve(2097152)
7.08 ms ± 80.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [561]: %timeit prime_wheel_sieve4(2097152)
7.1 ms ± 85.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [562]: %timeit prime_wheel_sieve4(4194304)
14.8 ms ± 120 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [563]: %timeit prime_wheel_sieve(4194304)
14.2 ms ± 145 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [564]: %timeit prime_wheel_sieve(8388608)
39.4 ms ± 1.44 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [565]: %timeit prime_wheel_sieve4(8388608)
41.7 ms ± 2.56 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
De acordo com meus testes, usar uma roda maior torna o código mais lento, e não mais rápido. Por que isso acontece? Teoricamente, usar uma roda maior estreita o espaço de busca, então não deveria aumentar o tempo de execução. Por que usar uma roda maior torna o código mais lento?
Ok, usando o método científico, controlei as diferenças entre as duas funções de modo que a única diferença entre elas, a única quantidade que pode afetar o desempenho, é a roda usada.
Fiz a primeira função usar o operador módulo em vez do AND bit a bit, e fiz a segunda função usar uma lista local, assim como a primeira. Eu queria separar código e dados, mas tudo bem.
@nb.njit(cache=True)
def prime_wheel_sieve3(n: int) -> np.ndarray:
wheel = [4, 2, 4, 2, 4, 6, 2, 6]
primes = np.ones(n + 1, dtype=np.uint8)
primes[:2] = False
for square, step in ((4, 2), (9, 6), (25, 10)):
primes[square::step] = False
k = 7
lim = int(math.sqrt(n) + 0.5)
i = 0
while k <= lim:
if primes[k]:
primes[k**2 :: 2 * k] = False
k += wheel[i]
i = (i + 1) % 8
return np.nonzero(primes)[0]
@nb.njit(cache=True)
def prime_wheel_sieve4_1(n: int) -> np.ndarray:
# fmt: off
wheel = [
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
]
# fmt: on
primes = np.ones(n + 1, dtype=np.uint8)
primes[:2] = False
for square, step in ((4, 2), (9, 6), (25, 10), (49, 14)):
primes[square::step] = False
k = 11
lim = int(math.sqrt(n) + 0.5)
i = 0
while k <= lim:
if primes[k]:
primes[k**2 :: 2 * k] = False
k += wheel[i]
i = (i + 1) % 48
return np.nonzero(primes)[0]
Tive que adicionar comentários de formatação para evitar que o Black Formatter bagunçasse minha tabela formatada no VS Code e, claro, isso não afeta o desempenho em nada.
As únicas diferenças entre as funções são o valor inicial de k, os primos que tiveram que ser processados antes de girar a roda propriamente dita e, claro, o tamanho da roda (e a própria roda). Estes tiveram que ser diferentes (porque usam rodas diferentes), mas todo o resto permanece inalterado.
In [679]: np.array_equal(prime_wheel_sieve3(65536), prime_wheel_sieve4_1(65536))
Out[679]: True
In [680]: %timeit prime_wheel_sieve3(65536)
162 μs ± 2.27 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [681]: %timeit prime_wheel_sieve4_1(65536)
158 μs ± 1.83 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [682]: %timeit prime_wheel_sieve3(131072)
326 μs ± 7.91 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [683]: %timeit prime_wheel_sieve4_1(131072)
322 μs ± 8.89 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [684]: %timeit prime_wheel_sieve3(262144)
659 μs ± 7.74 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [685]: %timeit prime_wheel_sieve4_1(262144)
655 μs ± 12.2 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [686]: %timeit prime_wheel_sieve3(524288)
1.45 ms ± 14.1 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [687]: %timeit prime_wheel_sieve4_1(524288)
1.42 ms ± 8.13 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In [688]: %timeit prime_wheel_sieve3(1048576)
3.2 ms ± 68.4 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [689]: %timeit prime_wheel_sieve4_1(1048576)
3.2 ms ± 199 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Curiosamente, prime_wheel_sieve4_1o desempenho é um pouco mais rápido que prime_wheel_sieve3, mas é só um pouquinho. A aceleração é insignificante, mas sei que o tempo de execução do código é estocástico por natureza, então prime_wheel_sieve4_1um desempenho consistentemente mais rápido que , prime_wheel_sieve3um pouquinho, é estatisticamente significativo. Embora eu não tenha testado muito, isso não exclui a possibilidade de coincidência.
Mas, de acordo com minha teoria, eu deveria ver uma redução de 14,29% no tempo de execução, e não basicamente nenhuma melhoria. Meus testes fortaleceram meu argumento.
Então por que usar uma roda maior não acelera o código?
O aumento de velocidade é insignificante porque você acelerou apenas uma parte insignificante do trabalho geral. A maior parte do tempo é gasta pelos
primes[...] = Falsecomandos, e eles são os mesmos para ambas as rodas: você faz isso para cada número primo até sqrt(n). Não importa que a roda maior encontre esses números primos um pouco mais rápido.