求出10的8次方以内的所有质数，看谁的代码最快,

[lizhongkun]

如上！

[Kaye]

这个问题我和LLnju在下面这个贴子上有详细讨论
http://www.csdn.net/expert/topic/417/417110.xml?temp=.3694422
intfree大虾的代码最快100亿不到10秒（统计个数无法知道详细的质数），大约60秒（知道所有质数）
LLnju的也比较快（具体时间我记不清楚了）
我的100亿大约需要200秒！！

[kbsoft]

那个是以6为基的素数表，算法较为复杂

[lizhongkun]

to: kaye非常佩服你，！！
你能解释一下intfree的pascal code吗？？
思想也行！！
我把它添出来：
{$R-,S-,Q-,O+}
{$APPTYPE CONSOLE}

program prime;

uses
sysutils;

const
maxlen3 = 1000;
dn = int64(10) shl 30;
dm = 7;

var
primes: array[1 .. maxlen3] of integer;
v: array of integer;
mark: array of boolean;
m, ans, Q, phiQ: integer;
n, len2, len3, sum: int64;
start: TDateTime;

procedure init;
var
l, k, s, max, sqmax: integer;
p: int64;
begin

max := round(exp(2 / 3 * ln(1.0 * n))) + 1;
sqmax := round(sqrt(max));
setlength(mark, max);
for s := 2 to sqmax do
if not mark[s] then
begin
k := s * s;
while k < max do
begin
mark[k] := true;
k := k + s;
end;
end;

l := 0; p := 2;
while p * p * p <= n do
begin
if not mark[p] then
begin
l := l + 1;
primes[l] := p;
end;
p := p + 1;
end;
len3 := l;

k := max - 1; sum := 0; s := 0;
while p * p <= n do
begin
if not mark[p] then
begin
l := l + 1;
while p * k > n do
begin
if not mark[k] then
s := s + 1;
k := k - 1;
end;
sum := sum + s;
end;
p := p + 1;
end;
len2 := l;

while p < max do
begin
if not mark[p] then
l := l + 1;
p := p + 1;
end;

sum := (len2 - len3) * l - sum;
sum := len3 * (len3 - 1) div 2 - len2 * (len2 - 1) div 2 + sum;

Q := 1;
phiQ := 1;
if len3 < dm then
m := len3
else
m := dm;
for s := 1 to m do
begin
Q := Q * primes[s];
phiQ := phiQ * (primes[s] - 1);
end;

setlength(v, Q);
for s := 0 to Q - 1 do
v[s] := s;
for s := 1 to m do
for k := Q - 1 downto 0 do
v[k] := v[k] - v[k div primes[s]];

end;

function phi(x: int64; a: integer): int64;
begin

if a = m then
begin
result := (x div Q) * phiQ + v[x mod Q];
exit;
end;

if x < primes[a] then
begin
result := 1;
exit;
end;

result := phi(x, a - 1) - phi(x div primes[a], a - 1);

end;

var
s: string;
begin
writeln('calc pi(n)');
write(format('n (default = %d) = ', [dn]));
try
readln(s);
n := strtoint64(s);
except
n := dn;
end;

start := now;
init;
ans := phi(n, len3) - sum + len3 - 1;

writeln(format('pi(%d) = %d (%.2fs)',
[n, ans, (now - start) * SecsPerDay]));
writeln('Press enter');
readln;
end.


请大家解释一下，！！
请大家解释一下思想！！！

等待更high的出现！！！！！
另，下面的ok 吗？？

#define MAX 100000000
long num[50000],*p,n=5;
int i,m=2;
num[0]=1;
num[1]=2;
num[2]=3;
p=num+3;
while(n<MAX)
{ for(i=2;i<=m;i++)
{ if(n%num[i]==0)
break;
else if(i==m) {*p=m;p++;m++}

}
n+=2;
}


明后天我会到机房调试，最快的奖励500分。

[NowCan]

#include <stdio.h>
#include <math.h>
#include <windows.h>

#undef _PRINT_RESULT_

void main()
{
unsigned long *check;
unsigned long num_check = 1;
unsigned long num_result = 2;
unsigned long number, counter;
unsigned long min, max, est, chkmax;
bool test;
min = 5;
printf("Input N:");
scanf("%lu",&max);
check=new unsigned long[sqrt(max)+1];
check[ 0 ] = 3;
#ifdef _PRINT_RESULT_
printf("2\t3\t");
#endif
est=GetTickCount();

for ( number = min; number <= max ; number+=2 )
{
test = true;
chkmax=sqrt( number );
for ( counter = 0; counter < num_check && *(check+counter) <= chkmax; ++counter )
{
if ( ( number % *(check+counter) ) == 0 )
{
test = false;
break;
}
}
if ( test )
{
#ifdef _PRINT_RESULT_
printf("%lu\t",number);
#endif
if ( number < sqrt( max ) )
*(check+num_check++) = number;
++num_result;
}
}
delete [] check;
printf("\n\nchecked:%lu Total:%lu\n",num_check,num_result);
printf("The program run for %lu ms\n",GetTickCount()-est);

}

这就是Kaye(菜菜狗) 的程序

[Kaye]

to nowcan:
这个程序经过你改编了吧，^_^！不过这不是我们最后讨论的结果，这个程序还比较慢，无法解答当初的100亿问题！我刚才找了一下，最后我的程序是下面这样的，虽然不怎么好，不过苟且能用了，呵呵
#include <iostream>
#include <cmath>
#include <windows.h>

using namespace std;

unsigned long check[10000];//用于存放小于10万的素数
unsigned long num_check = 1;

void initialize();
void prime();

void main()
{
int begin = ::GetTickCount();
initialize();
prime();
cout << "Using " << ((::GetTickCount() - begin)) << " ms" << endl;
}

//初始化100000以内的素数表
void initialize()
{
check[0] = 2;
const unsigned long min_check = 3;//2为偶数素数，所以从3开始
const unsigned long max_check = 100000;

for (unsigned long number = min_check; number <= max_check; number = number + 2)
{
unsigned long sqrt_num = sqrt(number);
bool test = true;
for (unsigned long i_check = 0; (i_check < num_check) && (check[i_check] <= sqrt_num); ++i_check)
{
if (number % check[i_check] == 0)
{
test = false;
break;
}
}
if (test)
{
check[num_check++] = number;
}
}
cout << num_check << endl;//素数表的总数
}

void prime()
{
const unsigned long min_seg = 1;//暂时还不能随意更改
const unsigned long max_seg = 1000;//与step的积应该大于10万(因为初始素数表里已经包含这些素数)
const unsigned long step = 100000;//按照理论10万大约有最好的效率,可以任意更改
//(测试发现50000是最快的，原因不明)
bool data[step + 1];//存储段的素数情况，需要和step大小一样
unsigned long pos_check[10000];//每段各素数最后的位置；
unsigned long total = 0;

total = num_check;
for (unsigned long i_pos = 0; i_pos < 10001; ++i_pos) pos_check[i_pos] = check[i_pos];

for (unsigned long i_seg = 1; i_seg <= max_seg; ++i_seg)
{
for (unsigned long i_check = 1; i_check < num_check; ++i_check)//此处可以从0开始(则输出total减一)
{ //则下面i_data+=2变成++i_data
unsigned long pos = pos_check[i_check];
unsigned long tmp_check = 2 * check[i_check];//乘2表示只检查素数的奇数倍
while (pos <= step)
{
data[pos] = false;
pos += tmp_check;
}
pos_check[i_check] = pos - step;
}

for (unsigned long i_data = 1; i_data <= step; i_data += 2)
{
if (data[i_data])
{
++total;
// cout << i_seg << '*' << step << '+' << i_data << endl;//打印实际素数的值
}
else
{
data[i_data] = true;//恢复data[i_data]初始值，以便下一轮计算
}

}
if (i_seg % 10 == 0)
cout << "Processing " << i_seg << '*' << step << ' ' << total << endl;
}
cout << total << endl;//若i_check从0开始，则将1 作为素数所以结果比实际多一,要减掉
}

[lizhongkun]

to :cui(蚊子王) ( ) 信誉：279 2002-05-19 11:50:00 得分:0


我打算用汇编编写一个，看看速度：）
我期待着，
上面的，我会好好测试的

[lbaby]

请各位考虑好你们的算法
算法对程序的影响可是致命的

[Kaye]

上面程序有一点错误
//cout << i_seg << '*' << step << '+' << i_data << endl;//打印实际素数的值
i_seg应该改为(i_seg - 1)
这个是输出实际素数值的错误（因为数据太多，这个输出实际上没有什么意义，当时也就没有特别检查，程序本身的计算应该是没有错的）

[lizhongkun]

知道了！

[kbsoft]

intfree高手的程序如下：
{$R-,S-,Q-,O+}
{$APPTYPE CONSOLE}

program prime;

const
maxn = int64(10) shl 30;
l = 1 shl 18;
ll = 10 shl (30 - 18);
maxlen = 10000;

var
primes, modprimes: array[1 .. maxlen] of integer;

procedure init;
var
i, k, p: integer;
begin
k := 1;
for i := 1 to maxlen do
repeat
k := k + 1; p := 1; primes := k;
while (primes[p] <= sqrt(k)) and (k mod primes[p] <> 0) do p
:= p + 1;
until (primes[p] = k) or (k mod primes[p] <> 0);
end;

var
sum, last: integer;
mark: array[0 .. l shr 1 - 1] of boolean;

procedure makeprimes;
var
i, k, s, p, sqmax: integer;
begin
sum := 0; fillchar(modprimes, sizeof(modprimes), $FF);
for i := 1 to ll do begin

fillchar(mark, sizeof(mark), true);
k := 1; sqmax := round(sqrt(1.0 * i * l));
repeat
k := k + 1; p := primes[k];
if p > sqmax then break;
if i > 1 then begin
if modprimes[k] = - 1 then
modprimes[k] := int64(i - 1) * l mod p
else
modprimes[k] := (modprimes[k] + l) mod p;
s := (p * (1 + modprimes[k] and 1) - modprimes[k]) shr 1;
end else s := p * p shr 1;
while s < l shr 1 do begin
mark[s] := false;
s := s + p;
end;
until false;
for k := 0 to l shr 1 - 1 do if mark[k] then sum := sum + 1;

if (int64(i) * l mod (10 shl 20) = 0) or (i = ll) then begin
writeln(int64(i) * l shr 20 : 5, 'M = ',
sum : 9, ' (', sum - last, ')');
last := sum;
end;

end;
end;

begin
init;
makeprimes;
end.

[kbsoft]

今天累了，回去给你写个C的

[kbsoft]

DEFLNG A-Z
DIM B(0 TO 15) AS INTEGER
DIM TM AS STRING
FOR I = 0 TO 15
READ B(I)
NEXT
DATA 1,2,4,8,16,32,64,128,256,512
DATA 1024,2048,4096,8192,16384,-32768
AAA:
INPUT "Input a number within 3-1000000:"; N
IF N < 3 OR N > 1536000 THEN
PRINT "Out of range!": GOTO AAA
END IF
TM = TIME$
K = INT(N / 48 + 1)
REDIM A(K) AS INTEGER
FOR I = 0 TO K
A(I) = 0
NEXT
A(0) = 1
Q = INT(SQR(N))
Q = INT((4 * Q - 3 - (Q MOD 6)) / 12)
H = INT((4 * N - 3 - (N MOD 6)) / 12)
FOR I = 1 TO Q
R = INT(I / 16): S = B(I MOD 16): T = A(R) AND S

IF T = 0 THEN

P = 3 * I + 1 + (I MOD 2)
M = INT((4 * P * P - 3 - (P * P MOD 6)) / 12)

U = 3 * I + 4 + ((I + 1) MOD 2)
V = INT((4 * P * U - 3 - ((P * U) MOD 6)) / 12)

W = V - M

FOR J = M TO H STEP 2 * P
R = INT(J / 16): S = B(J MOD 16): A(R) = A(R) OR S
R = INT((J + W) / 16)
IF R <= K THEN
S = B((J + W) MOD 16): A(R) = A(R) OR S
END IF
NEXT J
END IF
NEXT I
E = 2
PRINT USING "##########"; 2; 3;
FOR I = 1 TO H
R = INT(I / 16): S = B(I MOD 16): T = A(R) AND S
IF T = 0 THEN
P = 3 * I + 1 + (I MOD 2)
IF P > N THEN GOTO BBB
E = E + 1
PRINT USING "##########"; P;
END IF
NEXT I
BBB:
PRINT : PRINT "TOTAL:"; E
PRINT "START TIME:"; TM
PRINT "END TIME:"; TIME$
GOTO AAA
END

[wolaiye3]

我也考虑考虑！

[NowCan]

是啊，早就讨论过了，我还写了一个呢，只是速度不行。

[hyfyxu]

怎么没有人用C++编呀？？？
明天我把C++代码传上来。

[Kaye]

这个程序应该是60秒的程序吧，我对pascal也不太懂，所以这个程序也没有完全看懂，但是基本上都是基于那个贴子上的筛法实现的！如下：
1.一个数如果不能被小于其平方根的所有质数除那么它就是质数，如：11，小于sqrt(11)的质数有2，3。11不能被2，3除，那么11是质数
2.求小于n的所有质数，将小于sqrt(n)的质数列表的倍数全部除去，剩下的就是小于n的所有质数！如：小于100的所有质数，小于10的所有质数为2，3，5，7那么将小于100的2，3，5，7的倍数全部除去，剩下的就是小于100的质数
另外，当时考虑的是100亿，数据量太大，用了分段来解决这个问题！
intfree大虾的最后一个不到10s的程序，据他说用了容斥原理，我是完全不懂了^_^
顺便问一下 kbsoft(景乐)，你能解释一下你的第二个程序是什么意思吗?我实在是一点也没看懂！

[蚊子王]

我打算用汇编编写一个，看看速度：）

[Kaye]

to cui(蚊子王):
信誉：279
??????！！！！！！
你没有这么夸张吧？如何做到如此信誉的？

[zuoyefeng]

关注中。。。。。。。

[蚊子王]

还在讨论啊，考虑考虑