👨🏻‍🔬 ⚕️ ⛅️ 比Wikipedia更快筛选出数十亿个简单数字 👩🏼‍🤝‍👨🏾 👩🏻‍🏫 🉑

（ 图来源 ）

众所周知，Eratosthenes筛（RE）是最早出现在计算机发明之前的算法之一。因此，您可能会认为，经过数个世纪的研究，该算法已被反复研究，无法添加任何内容。如果您查看Wikipedia，则有大量关于权威资源的参考资料，您可以在其中轻松淹没。因此，当我无意中发现前一天在Wikipedia上显示为最佳选项时，我感到非常惊讶。

就是这样。在关于函数式编程（FP）的一篇文章的讨论中，他提出了一个问题：如何以这种范例编写RE。我不仅仅拥有很少的AF知识，我也不敢判断答案，但是讨论中的其他参与者拒绝了立即提出的一些解决方案，这表明，这不是理论上的复杂性

（ ）

$O（n \ log \ log n）$ （1）

拟议的实施将具有计算复杂性

（ ）

$O（n ^ 2 / \ log n）$ （2）

如此复杂，您迫不及待想要筛选出例如1000万个数字。我开始感兴趣，并尝试使用常规的程序编程来根据Wikipedia实现最佳版本。在Delphi-7中，我得到了以下代码：

清单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由具有相反值的筛布尔数组表示-如果数字为质数，则表示为false，这会减少筛选过程中否定运算（非）的次数。程序中包含3个过程：RE的初始化-初始化，计算（RE中数字的筛选和删除线）-calc和结果打印的质数的打印-打印，并对找到的数字进行计数。我将特别注意打印过程中注释掉的素数输出：对于在N = 1000的测试，注释将被删除。

在calc过程中，我使用Wikipedia建议：对于下一个素数i，从RE中删除数字

， ， ，

$i ^ 2，i ^ 2 + i，i ^ 2 + 2i，...$

该程序在17.6秒内筛选了十亿个数字。在我的PC（3.4 GHz Intel Core i7 CPU）上。
完成该程序后，我突然想起了偶数和奇数的众所周知的属性。

引理1。1 ）奇数+奇数=偶数； 2）奇数+偶数=奇数; 3）偶数+偶数=偶数。

证明

1）

（ ） （ ）

$（2n + 1）+（2m +1）= 2n + 2m + 2$ 可被2. TCD整除。
2）

（ ） （ ）

$（2n +1）+（2m）= 2n + 2m +1$ 不能被2整除。 chdd。
3）

（ ） （ ）

$（2n）+（2m）= 2n + 2m$ 可被2. TCD整除。

引理2。 一个奇数 的 平方是一个奇数。
证明。

（ ）

$（2n + 1）^ {2} = 4n ^ {2} + 4n + 1$ 不能被2整除。 chdd。

备注。 大于2的质数是奇数。

因此，您只能删除奇数：

$i ^ 2，i ^ 2 + 2i，i ^ 2 + 4i，...$ （3）

但是首先，您需要删除所有偶数。这是通过修改后的init初始化过程完成的。

清单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.

该程序运行了9.9秒。 -几乎快了一倍。

为了评估程序的实时操作与理论之间的对应关系，我建议

$dt = C * O，$

在哪里

$dt$ -测量的运行时间；

$C$ -与时间的维度一致；

$O$ -理论评估。

从这里计算

$C$ 评估（1）和（2）。对于

$N = 10 ^ 6$ 和更少

$dt$ 接近零。因此，我将数据带到第一个程序中

$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}$

如我们所见，估计（1）更接近实际结果。对于第二个程序，观察到类似的图片。我非常怀疑我是使用序列（3）发现美国的，我将非常感谢与使用此方法的工作的联系。 Wikipedia的作者很可能淹没在电子战的信息海中，但却错过了这项工作。

PS有关具有“线性运行时间”的Wikipedia算法，请参见

比Wikipedia更快筛选出数十亿个简单数字

More articles: