🔑 🧥 🌟 Tamizando miles de millones de números simples más rápido que Wikipedia 🏳️‍🌈 🍬 🔅

( Fuente de la figura )

Es bien sabido que el Tamiz de Eratóstenes (RE) es uno de los algoritmos más antiguos que aparecieron mucho antes de la invención de las computadoras. Por lo tanto, puede pensar que a lo largo de los siglos este algoritmo se ha estudiado de arriba abajo y no se le puede agregar nada. Si nos fijamos en Wikipedia, hay un mar de referencias a fuentes autorizadas en las que puede ahogarse fácilmente. Por lo tanto, me sorprendió cuando descubrí accidentalmente el otro día que la opción que se presenta como óptima en Wikipedia se puede optimizar significativamente.

Fue asi. En una discusión de un artículo sobre programación funcional (FP), hizo una pregunta : cómo escribir RE en este paradigma. Al poseer más que un escaso conocimiento de FP, no me atrevo a juzgar las respuestas, pero otros participantes en la discusión rechazaron algunas de las soluciones propuestas de inmediato, lo que indica que en lugar de la complejidad teórica

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

la implementación propuesta tendrá complejidad computacional

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

y que con tanta complejidad no puede esperar cuando, por ejemplo, se tamizan 10 millones de números. Me interesé e intenté implementar la versión óptima de acuerdo con Wikipedia, utilizando la programación de procedimientos habitual. En Delphi-7, obtuve el siguiente código:

Listado 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 está representado por la matriz booleana de tamiz con valores inversos: si el número es primo, se indica como falso, lo que reduce el número de operaciones de negación (no) durante el tamizado. Hay 3 procedimientos en el programa: inicialización del RE - init, cálculos (separación y tachado de números en el RE) - calc, e impresión de los números primos encontrados como resultado - print, mientras se cuenta el número de números encontrados. Especialmente llamaré la atención sobre la salida comentada de los números primos en el procedimiento de impresión: para probar en N = 1000, el comentario se elimina.

Aquí, en el procedimiento de cálculo, uso la recomendación de Wikipedia: para el próximo número primo i, elimine los números del RE

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

Este programa analizó mil millones de números en 17.6 segundos. en mi PC (CPU Intel Core i7 de 3.4 GHz).
Habiendo hecho este programa, de repente recordé las propiedades bien conocidas de los números pares e impares .

Lema 1. 1) impar + impar = par; 2) impar + par = impar; 3) par + par = par.

Prueba

1)

$(2n + 1) + (2m + 1) = 2n + 2m + 2$ divisible por 2. TCD.
2)

$(2n + 1) + (2m) = 2n + 2m + 1$ no divisible por 2 sin resto. Chtd.
3)

$(2n) + (2m) = 2n + 2m$ divisible por 2. TCD.

Lema 2. El cuadrado de un número impar es un número impar.
Prueba.

$(2n + 1) ^ {2} = 4n ^ {2} + 4n + 1$ no divisible por 2 sin resto. Chtd.

Observación Un número primo mayor que 2 es impar.

Por lo tanto, solo puede eliminar números impares:

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

Pero primero debes tachar todos los números pares. Esto se realiza mediante un procedimiento de inicialización init modificado.

Listado 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 funcionó durante 9.9 segundos. - Casi el doble de rápido.

Para evaluar la correspondencia de la operación en tiempo real del programa con la teoría, sugerí que

$dt = C * O,$

donde

$dt$ - tiempo de operación medido;

$C$ - constante con la dimensión del tiempo;

$O$ - evaluación teórica.

Calculado desde aquí

$C$ para evaluar (1) y (2). Para

$N = 10 ^ 6$ y menos

$dt$ cerca de cero Por lo tanto, traigo los datos del primer 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, la estimación (1) está mucho más cerca de los resultados reales. Para el segundo programa, se observa una imagen similar. Dudo mucho que descubrí América usando la secuencia (3) y estaré muy agradecido por el enlace al trabajo donde se aplicó este enfoque. Lo más probable es que los autores de Wikipedia se ahogaron en un mar de información sobre la guerra electrónica y se perdieron este trabajo.

PD Para el algoritmo de Wikipedia con "tiempo de ejecución lineal" ver

Tamizando miles de millones de números simples más rápido que Wikipedia

More articles: