庆五一，算法大比拼，放分大行动(300分)

使用类似压缩，查表 和 位操作，速度比我前面的稍快一点。基本思想同上一个是一样的。
const MMax = 20000000;
var
Bit: array [0..MMax div 30 + 1] of Byte


function CountPrimeNumber2: Integer;
const
Point: array [0..7] of Byte = (1,7,11,13,17,19,23,29);
Reflect: array [1..30] of Byte =
//1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
(0,0,0,0,0,0,1,1,1,1,2,2,3,3,3,3,4,4,5,5,5,5,6,6,6,6,6,6,7,7);
Offset: array [1..29] of Byte =
// 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
($01,0,0,0,0,0,$02,0,0,0,$04,0,$08,0,0,0,$10,0,$20,0,0,0,$40,0,0,0,0,0,$80);

var
ki, i, kj, j, vi, M, vv: Cardinal;
begin
FillChar(Bit, Sizeof(Bit), $00);
Bit[0] := $01
// 1 不处理
for ki := 0 to Round(Sqrt(MMax)/30) do
for i := 0 to 7 do
if (Bit[ki] and Offset[Point]) = 0 then
begin
vi := ki*30 + Point;
M := MMax div vi - 1;
kj := M div 30
//30 当成 29 处理
j := Reflect[M mod 30 + 1];
while (kj > ki) or ((kj = ki) and (j >= i)) do
begin
if (Bit[kj] and Offset[Point[j]]) = 0 then
begin
vv := vi * (kj * 30 + Point[j]);
Inc(Bit[vv div 30], Offset[vv mod 30]);
end;
if j > 0 then
Dec(j)
else begin
j := 7;
Dec(kj);
end;
end;
end;
Result := 3
//2, 3, 5
for ki := 0 to MMax div 30 - 1 do
for i := 0 to 7 do
if (Bit[ki] and Offset[Point]) = 0 then //Num = ki * 30 + Point
Inc(Result);
ki := MMax div 30;
for i := 0 to 7 do
begin
if ki * 30 + Point > MMax then Exit;
if (Bit[ki] and Offset[Point]) = 0 then //Num = ki * 30 + Point
Inc(Result);
end;
end;

不用压缩位的话，可以直接用 数组操作，速度差不多。
function CountPrimeNumber3: Integer;
const
Point: array [0..7] of Byte = (1, 7,11,13,17,19,23,29);
Step: array [0..7] of Byte = (6, 4, 2, 4, 2, 4, 6, 2);
Reflect: array [1..30] of Byte =
//1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
(0,0,0,0,0,0,1,1,1,1,2,2,3,3,3,3,4,4,5,5,5,5,6,6,6,6,6,6,7,7);

var
ki, i, j, vi, vj, Mo: Cardinal;
begin
FillChar(B, Sizeof(B), True);
B[1] := False
// 1 不处理
for ki := 0 to Round(Sqrt(MMax)/30) do
for i := 0 to 7 do
begin
vi := ki*30 + Point;
if not B[vi] then Continue;

vj := MMax div vi;
Mo := (vj - 1) mod 30 + 1;
j := Reflect[Mo];
vj := vj + Point[j] - Mo;
while vj >= vi do
begin
if B[vj] then B[vi * vj] := False;
if j > 0 then Dec(j) else j := 7;
Dec(vj, Step[j]);
end;
end;
Result := 3
//2, 3, 5
for ki := 0 to MMax div 30 - 1 do
for i := 0 to 7 do
if B[ki * 30 + Point] then Inc(Result);
ki := MMax div 30;
for i := 0 to 7 do
begin
if ki * 30 + Point > MMax then Exit;
if B[ki * 30 + Point] then Inc(Result);
end;
end;

我已经回到家了
感谢大家的参与和讨论
希望以下的探讨围绕以下两方面
1.算法的改进追求最快的速度
2.能够提出更多的实现方法，这能扩展我们的视野。

希望这个帖子可以给我和大家在编程深度和广度上都有帮助
祝大家五一快乐

路过的朋友记得要收藏这个帖子啊：）

看到下面的话，我来回上几句，说说Prime的算法原理，Prime也是目前已知用除法最快的算法了：
===============================================================================
来自：DoubleWood, 时间：2006-5-3 0:53:40, ID：3433017
......
目前见到最快的应该是PrimeNumber.zip里的算法，但是没有源码，无法得知算法。
Delphi写的我见到最快的是Prime.dcu，作者同样没有提供源码，但是提供了声明部分，所以可可以在程序里调用，大约比jeffrey_s的算法快上几十倍。作者是Hagen Reddmann。著名的DEC（ Delphi Encryption Compendium ）的作者。
汇编肯定要比高级语言要快，前提是你的代码正确[]，但是提升未必有想象中的高。很多时候，一个好的算法比用汇编重写的提升要来得快
===============================================================================
先说说一般的求质数算法实现：
就是依次循环变量I从2到N，然后检查它是不是质数，如果是就显示，不是就检查下一个。
问题是这个方法会做很多重复的检查，所以程序的速度会变慢。
比如2是质数，但是2的倍数都不是质数，如果从让I从2变到N，就会多测试4/6/8/10...
所以I：最好只测试2和所有的奇数就可以了。
II：同理，可以把所有3的倍数也排除在测试之内
综合I和II，可以的出，在任意连续的6个数[6J / 6J+1 / 6J+2 / 6J+3 / 6J+4 / 6J+5]中，就只用测试两个而已[6J+1 / 6J+5]，为什么？自己想去，因此算法工作量就是原来算法的2/6 = 1/3倍。

再往下走，假设I是上面的一个数[非2/3的倍数，即：I = 6J+1 / 6J+5]，如何测试I是质数呢？用的着依次用2/3/4.....I来除测试看能不能整除吗？
首先来说， 不用，只需要测试2/3/4.....(I/2)就可以，道理很简单
在分析下去，还可以缩小范围成2/3/4.....Sqrt(I)<对I开平方>

但是，这个也是浪费，最好的方法就是用质数去除I，因为之前已经除去了2，3的倍数，所以就从5到 Sqrt(I)的质数来除I即可。

大体算法就是这样了,给出如下代码,基于C的， 很简单，核心代码超不过十行，自己转换吧:

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
/* ------------------------------------------------------ */
/* PROGRAM prime number generator : */
/* This program generates first MAXSIZE-1 primes by */
/* using the conventional division method, but multiples */
/* of 2 and 3 are not tested. Therefore it is faster. */
/* */
/* Copyright Ching-Kuang Shene July/02/1989 */
/* ------------------------------------------------------ */
#include <stdio.h>
#define MAXSIZE 100
#define YES 1
#define NO 0
void main(void)
{
int prime[MAXSIZE]
/* array to contains primes */
int gap = 2
/* increasing gap = 2 and 4 */
int count = 3
/* no. of primes */
int may_be_prime = 5
/* working variable */
int i, is_prime;

prime[0] = 2
/* Note that 2, 3 and 5 are */
prime[1] = 3
/* primes. */
prime[2] = 5;
while (count < MAXSIZE)
{
/* loop until prime[] full*/
may_be_prime += gap
/* get next number */
gap = 6 - gap
/* switch to next gap */
is_prime = YES
/* suppose it were a prime*/
for (i = 2
prime*prime <= may_be_prime && is_prime
i++)
{
if (may_be_prime % prime == 0) /* NO */
is_prime = NO
/* exit */
}
if (is_prime) /* if it IS a prime... */
prime[count++] = may_be_prime
/* save it */
}

for (i = 0
i < count
i++)
{
if (i % 10 == 0) printf(&quot;/n&quot;
printf(&quot;%5d&quot;, prime);
}
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

如上给出的算法应该是给出的最快的了吧^-^

现在不在家,没有环境啊
楼主,帮你顶了

很遗憾，网事如风给出的算法非常的慢……循环次数可能少了，但是每次循环里面都作了
一个非常耗时的操作 mod.而前面大多数算法循环里大部分做的仅是判断而已，或者是乘
法.而且你说的算法根本就不是Prime.dcu里的算法。
又：
jeffrey_s最后给出的算法还有改进余地，不过提升已经很少了，
vi := ki*30 + Point;
if B[ki * 30 + Point] then Inc(Result);
这两处都可以改成加法，在我机器上大约提升2~3%左右吧。

唉,你们慢慢算吧,我没空

翻译了一下 DoubleWood 所介绍的算法。

var
BBase: array [0..$FFFF] of Boolean;
BLine: array [0..$FFFF] of Boolean;
BBackup: array [0..$FFFF] of Boolean;

function CountPrimeNumber5: Integer;
const
CACHE: Cardinal = $10000
//65536;
var
i, j, k, v, Qrt: Cardinal;
begin
for i := 0 to CACHE - 1 do
BBackup := Odd(i);

Qrt := Round(Sqrt(MMax));
Move(BBackup, BBase, Qrt);
BBase[2] := True;
for i := 3 to Round(Sqrt(Qrt)) do //Use 1..Sqrt(Qrt) to Find 1..Qrt
if BBase then
begin
j := i * i
//for ji := i to Qrt div i do
repeat
BBase[j] := False
//B[i * ji] := False;
Inc(j, i);
until (j > Qrt);
end;

Result := 0;
for k := 0 to MMax div CACHE do
begin
Move(BBackup, BLine, CACHE);
v := k * CACHE;
for i := 3 to Round(Sqrt(v + CACHE)) do
if BBase then
begin
if k = 0 then //if k = 0 save Base
j := i * i
else //else Let in [v..v+Cache];
j := (v div i + 1) * i - v;
while j < CACHE do
begin
BLine[j] := False
//B[i * ji] := False;
Inc(j, i);
end;
end;
for i := 0 to Length(BLine) - 1 do
begin
if v + i > MMax then Exit;
if BLine then Inc(Result);
end;
end;
end;

压缩偶数位后变为：
function CountPrimeNumber6: Integer;
const
CACHE = $8000;
var
ODDBase, ODDLine: array [0..CACHE-1] of Boolean;
var
i, j, PrimeNum, k, v, Qrt, M: Cardinal;
begin
Qrt := Round(Sqrt(MMax)) div 2;
FillChar(ODDBase, Qrt, True);
for i := 1 to Round(Sqrt(Qrt/2)) do
if ODDBase then
begin
PrimeNum := i * 2 + 1;
j := (PrimeNum + 1) * i
// = PrimeNum * PrimeNum div 2;
repeat
ODDBase[j] := False;
Inc(j, PrimeNum);
until (j > Qrt);
end;

Result := 1
// 2
M := MMax div 2;
for k := 0 to M div CACHE do
begin
FillChar(ODDLine, CACHE, True);
if k = 0 then ODDLine[0] := False
//1 is not a prime Num
v := k * CACHE;
for i := 1 to Round(Sqrt((v + CACHE)*2)) do
if ODDBase then
begin
PrimeNum := i * 2 + 1;
if k = 0 then
j := (PrimeNum + 1) * i // = PrimeNum * PrimeNum div 2
else
j := ((v * 2 div PrimeNum + 1) div 2 * 2 + 1) * PrimeNum div 2 - v;
while j < CACHE do
begin
ODDLine[j] := False;
Inc(j, PrimeNum);
end;
end;
for i := 0 to CACHE - 1 do
begin
if v + i > M then Break;
if ODDLine then Inc(Result);
end;
end;
end;

上上面的 j := (v div i + 1) * i - v;
改为 ((v div i + 1) div 2 * 2 + 1) * i - v;
那么原来 Inc(j, i)
就能改成 Inc(j, i * 2);类似的。减少了偶数的对比次数。

晕倒，还不结帖~！我取消订阅了

结帖~！

我的10000000要112.109秒的吗。
包含把算出的质数放到LISTBOX中。

楼上的把代码贴出来啊

很多算法都很优秀
看来我要食言了
基本上平均分配积分吧
重在学习