🙇🏻 🌜 🤽🏽 Selecionando bilhões de números simples mais rapidamente que a Wikipedia ⬆️ 🦀 👨🏽‍✈️

( Fonte da Figura )

É sabido que a Peneira de Eratóstenes (ER) é um dos algoritmos mais antigos que surgiram muito antes da invenção dos computadores. Portanto, você pode pensar que, ao longo dos séculos, esse algoritmo foi estudado de cima para baixo e nada pode ser adicionado a ele. Se você olhar para a Wikipedia, há um mar de referências a fontes autorizadas nas quais você pode se afogar facilmente. Portanto, fiquei surpreso quando descobri acidentalmente outro dia que a opção apresentada como ideal na Wikipedia pode ser significativamente otimizada.

Foi assim. Em uma discussão de um artigo sobre programação funcional (FP), ele fez uma pergunta : como escrever RE neste paradigma. Possuindo um conhecimento mais do que escasso da FA, não me atrevo a julgar as respostas, mas outros participantes da discussão rejeitaram algumas das soluções propostas imediatamente, indicando que, em vez de complexidade teórica

$O (n \ log \ log n)$ (1)

a implementação proposta terá complexidade computacional

$O (n ^ 2 / \ log n)$ 2)

e que com tanta complexidade você não pode esperar quando, por exemplo, 10 milhões de números são analisados. Fiquei interessado e tentei implementar a versão ideal de acordo com a Wikipedia, usando a programação processual usual. No Delphi-7, obtive o seguinte código:

Listagem 1

program EratosthenesSieve; // Sieve of Eratosthenes {$APPTYPE CONSOLE} uses SysUtils, DateUtils,Math; const version ='1.0.1d1'; N = 1000000000; // number of primes var sieve : array [2..N] of boolean; // false if prime t0, t1,dt : TDateTime; O,C : Extended; procedure init; var i : integer; begin for i:=2 to n do sieve [i] := false; end; //init procedure calc (start : integer); var prime, i : integer; breakLoop, exitProc : Boolean; begin prime := start; exitProc := false; repeat // find next prime prime := prime+1; while (prime<N) and sieve[prime] do inc (prime); i:= sqr(prime); // delete if i<= N then begin breakLoop := false; repeat if i<= N then begin sieve [i] := true; inc (i,prime); end else // if i<= N breakLoop := true; until breakLoop; end else // if prime+prime<= N exitProc := true; until exitProc; end; //calc procedure print; var i :integer; found : integer; begin found := 0; for i:=2 to N do if not sieve [i] then begin // write (i,', '); inc(found); end; writeln; writeln ('Found ',found,' primes.'); end; // begin // program body writeln ('Sieve of Eratosthenes version ', version); writeln('N= ',N); init; t0 := now; writeln('Program started ',DateTimeToStr(t0)); t0 := now; calc (1); t1 := now; writeln('Program finished ',DateTimeToStr(t1)); dt := SecondSpan(t1,t0); writeln ('Time is ',FloatToStr(dt),' sec.'); O := N* ln(ln(N)); C := dt/O; writeln ('O(N ln ln N)= ',O,' C=',C); O := sqr(N)/ln(N); C := dt/O; writeln ('O(N*N/ln N)= ',O,' C=',C); print; writeln ('Press Enter to stop...'); readln; end.

RE é representado pela matriz booleana da peneira com valores inversos - se o número for primo, é indicado como falso, o que reduz o número de operações de negação (não) durante a peneiração. Existem três procedimentos no programa: inicialização do RE - init, cálculos (peneiração e tacada de números no RE) - calc e impressão dos números primos encontrados como resultado - print e o número de números encontrados é contado. Vou chamar especialmente a atenção para a saída comentada de números primos no procedimento de impressão: para testes em N = 1000, o comentário é removido.

Aqui no procedimento de cálculo, uso a recomendação da Wikipedia: para o próximo número primo i, exclua os números do RE

$i ^ 2, i ^ 2 + i, i ^ 2 + 2i, ...$

Esse programa analisou um bilhão de números em 17,6 segundos. no meu PC (CPU Intel Core i7 de 3,4 GHz).
Tendo criado este programa, de repente me lembrei das propriedades conhecidas dos números pares e ímpares .

Lema 1. 1) ímpar + ímpar = par; 2) ímpar + par = ímpar; 3) par + par = par.

Prova

1)

$(2n + 1) + (2m + 1) = 2n + 2m + 2$ divisível por 2. TCD.
2)

$(2n + 1) + (2m) = 2n + 2m + 1$ não divisível por 2 sem restante. Chtd.
3)

$(2n) + (2m) = 2n + 2m$ divisível por 2. TCD.

Lema 2. O quadrado de um número ímpar é um número ímpar.
Prova.

$(2n + 1) ^ {2} = 4n ^ {2} + 4n + 1$ não divisível por 2 sem restante. Chtd.

Observação. Um número primo maior que 2 é ímpar.

Portanto, você pode excluir apenas números ímpares:

$i ^ 2, i ^ 2 + 2i, i ^ 2 + 4i, ...$ (3)

Mas primeiro você precisa cruzar todos os números pares. Isso é feito por um procedimento de inicialização init modificado.

Listagem 2

 program EratosthenesSieve; // Sieve of Eratosthenes {$APPTYPE CONSOLE} uses SysUtils, DateUtils,Math; const version ='1.0.1d1'; N = 1000000000; // number of primes var sieve : array [2..N] of boolean; // false if prime t0, t1,dt : TDateTime; O,C : Extended; procedure init; var i : integer; begin for i:=2 to n do sieve [i] := not odd(i); end; //init procedure calc (start : integer); var prime,prime2, i : integer; breakLoop, exitProc : Boolean; begin prime := start; exitProc := false; repeat // find next prime prime := prime+1; while (prime<N) and sieve[prime] do inc (prime); // i:= prime*prime; i:= sqr(prime); prime2 := prime+prime; // delete if i<= N then begin breakLoop := false; repeat if i<= N then begin sieve [i] := true; inc (i,prime2); end else // if i<= N breakLoop := true; until breakLoop; end else // if prime+prime<= N exitProc := true; until exitProc; sieve [2] := false; end; //calc procedure print; var i :integer; found : integer; begin found := 0; for i:=2 to N do if not sieve [i] then begin // write (i,', '); inc(found); end; writeln; writeln ('Found ',found,' primes.'); end; // begin // program body writeln ('Sieve of Eratosthenes version ', version); writeln('N= ',N); init; t0 := now; writeln('Program started ',DateTimeToStr(t0)); t0 := now; calc (2); t1 := now; writeln('Program finished ',DateTimeToStr(t1)); dt := SecondSpan(t1,t0); writeln ('Time is ',FloatToStr(dt),' sec.'); O := N* ln(ln(N)); C := dt/O; writeln ('O(N ln ln N)= ',O,' C=',C); O := sqr(N)/ln(N); C := dt/O; writeln ('O(N*N/ln N)= ',O,' C=',C); print; writeln ('Press Enter to stop...'); readln; end.

Este programa funcionou por 9,9 segundos. - quase o dobro da velocidade.

Para avaliar a correspondência da operação em tempo real do programa com a teórica, sugeri que

$dt = C * O,$

onde

$dt$ - tempo de operação medido;

$C$ - constante com a dimensão do tempo;

$O$ - avaliação teórica.

Calculado a partir daqui

$C$ avaliar (1) e (2). Para

$N = 10 ^ 6$ e menos

$dt$ perto de zero. Por isso, trago os dados do primeiro programa para grandes

$N.$

$N$	(1)	2)
$10 ^ 7$	$1,69 \ cdot 10 ^ {- 9}$	$2,74 \ cdot 10 ^ {- 9}$
$10 ^ 8$	$5,14 \ cdot 10 ^ {- 9}$	$1,47 \ cdot 10 ^ {- 8}$
$10 ^ 9$	$5,80 \ cdot 10 ^ {- 9}$	$1,29 \ cdot 10 ^ {- 7}$

Como vemos, a estimativa (1) está muito mais próxima dos resultados reais. Para o segundo programa, uma imagem semelhante é observada. Duvido muito que tenha descoberto a América usando a sequência (3) e ficarei muito grato pelo link para o trabalho em que essa abordagem foi aplicada. Provavelmente, os próprios autores da Wikipedia se afogaram em um mar de informações sobre a guerra eletrônica e perderam esse trabalho.

PS Para o algoritmo da Wikipedia com "tempo de execução linear", consulte

Selecionando bilhões de números simples mais rapidamente que a Wikipedia

More articles: