🤽🏾 ⭕️ 🚇 Filtrer des milliards de nombres simples plus rapidement que Wikipedia 👨‍🍳 🐻 👨🏽‍🤝‍👨🏻

( Source de la figure )

Il est bien connu que le tamis d'Ératosthène (RE) est l'un des algorithmes les plus anciens apparus bien avant l'invention des ordinateurs. Par conséquent, vous pourriez penser qu'au cours des siècles, cet algorithme a été étudié de haut en bas et que rien ne peut y être ajouté. Si vous regardez Wikipedia, il y a une mer de références à des sources faisant autorité dans lesquelles vous pouvez facilement vous noyer. Par conséquent, j'ai été surpris lorsque j'ai découvert accidentellement l'autre jour que l' option présentée comme optimale sur Wikipédia pouvait être considérablement optimisée.

C'était comme ça. Dans une discussion d'un article sur la programmation fonctionnelle (FP), il a posé une question : comment écrire RE dans ce paradigme. Possédant plus qu'une maigre connaissance de la FA, je n'ose pas juger les réponses, mais d'autres participants à la discussion ont rejeté immédiatement certaines des solutions proposées, indiquant qu'au lieu de la complexité théorique

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

la mise en œuvre proposée aura une complexité de calcul

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

et qu'avec une telle complexité, vous ne pouvez pas attendre lorsque, par exemple, 10 millions de numéros sont triés. Je me suis intéressé et j'ai essayé d'implémenter la version optimale selon Wikipedia, en utilisant la programmation procédurale habituelle. Dans Delphi-7, j'ai obtenu le code suivant:

Listing 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 représenté par le tableau booléen tamis avec des valeurs inverses - si le nombre est premier, il est indiqué comme faux, ce qui réduit le nombre d'opérations de négation (pas) pendant le tamisage. Il y a 3 procédures dans le programme: initialisation du RE - init, calculs (criblage et barrage des nombres dans le RE) - calc, et impression des nombres premiers trouvés en conséquence - print, et le nombre de nombres trouvés est compté. J'attirerai particulièrement l'attention sur la sortie commentée des nombres premiers dans la procédure d'impression: pour les tests à N = 1000, le commentaire est supprimé.

Ici, dans la procédure de calcul, j'utilise la recommandation Wikipedia: pour le prochain nombre premier i, supprimez les nombres de l'ER

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

Ce programme a passé au crible un milliard de chiffres en 17,6 secondes. sur mon PC (processeur Intel Core i7 à 3,4 GHz).
Ayant réalisé ce programme, je me suis soudain souvenu des propriétés bien connues des nombres pairs et impairs .

Lemme 1. 1) impair + impair = pair; 2) impair + pair = impair; 3) pair + pair = pair.

Preuve

1)

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

$(2n + 1) + (2m) = 2n + 2m + 1$ non divisible par 2 sans reste. Chtd.
3)

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

Lemme 2. Le carré d'un nombre impair est un nombre impair.
Preuve.

$(2n + 1) ^ {2} = 4n ^ {2} + 4n + 1$ non divisible par 2 sans reste. Chtd.

Remarque. Un nombre premier supérieur à 2 est impair.

Par conséquent, vous ne pouvez supprimer que des nombres impairs:

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

Mais vous devez d'abord barrer tous les nombres pairs. Cela se fait par une procédure d'initialisation init modifiée.

Listing 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.

Ce programme a fonctionné pendant 9,9 secondes. - presque deux fois plus vite.

Pour évaluer la correspondance entre le fonctionnement du programme en temps réel et le fonctionnement théorique, j'ai suggéré que

$dt = C * O,$

où

$dt$ - durée de fonctionnement mesurée;

$C$ - constant avec la dimension du temps;

$O$ - évaluation théorique.

Calculé à partir d'ici

$C$ évaluer (1) et (2). Pour

$N = 10 ^ 6$ et moins

$dt$ proche de zéro. Par conséquent, j'apporte les données sur le premier programme pour les grands

$N.$

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

Comme nous le voyons, l'estimation (1) est beaucoup plus proche des résultats réels. Pour le deuxième programme, une image similaire est observée. Je doute beaucoup d'avoir découvert l'Amérique en utilisant la séquence (3) et je serai très reconnaissant du lien avec le travail où cette approche a été appliquée. Très probablement, les auteurs de Wikipédia eux-mêmes se sont noyés dans une mer d'informations sur la guerre électronique et ont raté ce travail.

PS Pour l'algorithme Wikipedia avec «runtime linéaire», voir

Filtrer des milliards de nombres simples plus rapidement que Wikipedia

More articles: