33,008
社区成员
发帖
与我相关
我的任务
分享擂台：筛法求质数Sieve of Eratosthenes
筛法求质数Sieve of Eratosthenes描述：
http://primes.utm.edu/glossary/page.php?sort=SieveOfEratosthenes
相关帖子:
擂台：求N!的最后九个非零尾数
http://expert.csdn.net/Expert/topic/2986/2986867.xml?temp=.511532
＃＃谁愿意帮我做个测试，最好是P4的CPU！＃＃
http://expert.csdn.net/Expert/topic/3027/3027142.xml?temp=.850918
程序下载（最新版，暂不提供源码，仍在不断改进中）：
http://www.77studio.net/files/PrimeLabPro.zip
参考源代码：
http://www.77studio.net/files/PrimeLab_Src.zip
http://www.77studio.net/files/PrimeLab_Src_New.zip （比前一个新）
另一种更快的方法：
Sieve of Atkins
介绍：http://primes.utm.edu/links/programs/sieves/binary_quadratic/index.html
程序：http://www.77studio.net/files/SieveAtkins.zip
特别感谢xstring(麻雀), liangbch(宝宝), languagec(各有所求)（是从你的帖子开始研究质数的^_^）等，只为交流，适当的时候会和大家分享代码。
注：
以上未嵌入汇编，未使用ICC(Intel C++ Compiler)优化，VS 6.0 + SP5编译
http://primes.utm.edu/glossary/page.php?sort=SieveOfEratosthenes
相关帖子:
擂台：求N!的最后九个非零尾数
http://expert.csdn.net/Expert/topic/2986/2986867.xml?temp=.511532
＃＃谁愿意帮我做个测试，最好是P4的CPU！＃＃
http://expert.csdn.net/Expert/topic/3027/3027142.xml?temp=.850918
程序下载（最新版，暂不提供源码，仍在不断改进中）：
http://www.77studio.net/files/PrimeLabPro.zip
参考源代码：
http://www.77studio.net/files/PrimeLab_Src.zip
http://www.77studio.net/files/PrimeLab_Src_New.zip （比前一个新）
另一种更快的方法：
Sieve of Atkins
介绍：http://primes.utm.edu/links/programs/sieves/binary_quadratic/index.html
程序：http://www.77studio.net/files/SieveAtkins.zip
特别感谢xstring(麻雀), liangbch(宝宝), languagec(各有所求)（是从你的帖子开始研究质数的^_^）等，只为交流，适当的时候会和大家分享代码。
注：
以上未嵌入汇编，未使用ICC(Intel C++ Compiler)优化，VS 6.0 + SP5编译
...全文
请发表友善的回复…
发表回复
shines77 2004-09-14
- 打赏
- 举报
按照CSDN的精神，此贴最近将结
liangbch 2004-08-12
- 打赏
- 举报
我改进了我的程序，但仍然使用字节数组，分段筛，每段为 262140，筛4*10^9以内的所有素数，筛过程用时15.38秒，统计用时4.26秒，共计20.10秒(有一部分计算时间不在前2者之内).
若改为每段120K,总时间为22.60秒.
楼主的程序SieveAtkins 为在我的电脑上运行速度为6.32秒。
字节数组最大的问题在于读写内存的数据量较大，特别是统计素数的个数部分，很难提高速度。
若改为每段120K,总时间为22.60秒.
楼主的程序SieveAtkins 为在我的电脑上运行速度为6.32秒。
字节数组最大的问题在于读写内存的数据量较大，特别是统计素数的个数部分，很难提高速度。
NowCan 2004-08-10
- 打赏
- 举报
如果秋初素数只要50秒的话，那求个数绝对不用10秒。
DiabloWalkOnTheEarth 2004-08-10
- 打赏
- 举报
好像intfree的程序求10G内素数在 C2 535上求个数只要<10秒,求出素数只要50秒吧,咔咔
liangbch 2004-08-10
- 打赏
- 举报
我使用字节数组，分段筛，每段为 65520 ，筛4*10^9以内的所有素数，筛过程用时22秒，统计用时6秒，共计28.5秒。楼主的程序SieveAtkins 为6.32秒。
DiabloWalkOnTheEarth 2004-08-10
- 打赏
- 举报
2: 求素数个数的,咔咔,好像以前的大虾如 starfish , intfree, mathe 都不来了啊
{$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.
{$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.
DiabloWalkOnTheEarth 2004-08-10
- 打赏
- 举报
这是 intfree 的程序 :
1 : 可以求出具体素数的 ;
{$APPTYPE CONSOLE}
{$R-,S-,Q-,O-}
program prime;
uses sysutils, windows;
const
maxn = int64(10) shl 30;
maxlen = 10000;
m = 6;
dots = 80;
type
Pls = ^Tls;
Tls = array[byte] of byte;
var
mask: pointer;
primes, mprimes, ans: array[1 .. maxlen] of integer;
bits: array[byte] of integer;
len, r: integer;
start: TDateTime;
procedure solve(const a, b, c: integer; var x, y: integer);
begin
if a = 0 then begin
x := 0; y := c div b;
end else begin
solve(b mod a, a, c, y, x);
x := x - b div a * y;
end;
end;
procedure init;
var
p, k: integer;
begin
k := 2; len := 0;
repeat
len := len + 1; primes[len] := k;
repeat
k := k + 1; p := 1;
while (sqr(primes[p]) < k) and (k mod primes[p] <> 0) do p := p + 1;
until k mod primes[p] <> 0;
until 1.0 * k * k > maxn;
bits[0] := 0; for k := 1 to $FF do bits[k] := bits[k shr 1] + k and 1;
end;
procedure makeprimes;
var
i, k, n, size, kmax, min, last, p, x, dot: integer;
begin
n := 1; for i := 1 to m do n := n * primes[i];
for i := m + 1 to len do begin
p := primes[i]; solve(p, - n, 1, x, mprimes[i]);
mprimes[i] := (mprimes[i] mod p + p) mod p;
end;
r := len; size := (maxn div n + 7) shr 3; last := 0; dot := 0;
fillchar(ans, sizeof(ans), 0); getmem(mask, size);
for i := 1 to n - 1 do begin
k := 1; while (k < m) and (i mod primes[k] <> 0) do k := k + 1;
if i mod primes[k] = 0 then continue;
p := i - last; last := i; kmax := (maxn - i) div n;
for k := m + 1 to len do ans[k] := (ans[k] + mprimes[k] * p) mod primes[k];
fillchar(mask^, size, 0);
for k := m + 1 to len do begin
p := primes[k];
if ans[k] = 0 then min := p - 1 else min := ans[k] - 1;
asm
mov EAX, mask
mov EBX, kmax
mov ECX, min
mov EDX, p
@001:
BTS [EAX], ECX
add ECX, EDX
cmp ECX, EBX
jl @001
end;
end;
for k := 0 to size - 1 do kmax := kmax - bits[Pls(mask)^[k]];
r := r + kmax;
if i / n > dot / dots then begin
dot := dot + 1; write('.');
end;
end;
freemem(mask, size);
writeln;
end;
begin
SetPriorityClass(GetCurrentProcess, REALTIME_PRIORITY_CLASS);
start := now;
init;
makeprimes;
messageboxA(0,
PChar(Format('sum = %d, Time = %.2f',
[r, (now - start) * SecsPerDay])),
'prime', 0);
end.
1 : 可以求出具体素数的 ;
{$APPTYPE CONSOLE}
{$R-,S-,Q-,O-}
program prime;
uses sysutils, windows;
const
maxn = int64(10) shl 30;
maxlen = 10000;
m = 6;
dots = 80;
type
Pls = ^Tls;
Tls = array[byte] of byte;
var
mask: pointer;
primes, mprimes, ans: array[1 .. maxlen] of integer;
bits: array[byte] of integer;
len, r: integer;
start: TDateTime;
procedure solve(const a, b, c: integer; var x, y: integer);
begin
if a = 0 then begin
x := 0; y := c div b;
end else begin
solve(b mod a, a, c, y, x);
x := x - b div a * y;
end;
end;
procedure init;
var
p, k: integer;
begin
k := 2; len := 0;
repeat
len := len + 1; primes[len] := k;
repeat
k := k + 1; p := 1;
while (sqr(primes[p]) < k) and (k mod primes[p] <> 0) do p := p + 1;
until k mod primes[p] <> 0;
until 1.0 * k * k > maxn;
bits[0] := 0; for k := 1 to $FF do bits[k] := bits[k shr 1] + k and 1;
end;
procedure makeprimes;
var
i, k, n, size, kmax, min, last, p, x, dot: integer;
begin
n := 1; for i := 1 to m do n := n * primes[i];
for i := m + 1 to len do begin
p := primes[i]; solve(p, - n, 1, x, mprimes[i]);
mprimes[i] := (mprimes[i] mod p + p) mod p;
end;
r := len; size := (maxn div n + 7) shr 3; last := 0; dot := 0;
fillchar(ans, sizeof(ans), 0); getmem(mask, size);
for i := 1 to n - 1 do begin
k := 1; while (k < m) and (i mod primes[k] <> 0) do k := k + 1;
if i mod primes[k] = 0 then continue;
p := i - last; last := i; kmax := (maxn - i) div n;
for k := m + 1 to len do ans[k] := (ans[k] + mprimes[k] * p) mod primes[k];
fillchar(mask^, size, 0);
for k := m + 1 to len do begin
p := primes[k];
if ans[k] = 0 then min := p - 1 else min := ans[k] - 1;
asm
mov EAX, mask
mov EBX, kmax
mov ECX, min
mov EDX, p
@001:
BTS [EAX], ECX
add ECX, EDX
cmp ECX, EBX
jl @001
end;
end;
for k := 0 to size - 1 do kmax := kmax - bits[Pls(mask)^[k]];
r := r + kmax;
if i / n > dot / dots then begin
dot := dot + 1; write('.');
end;
end;
freemem(mask, size);
writeln;
end;
begin
SetPriorityClass(GetCurrentProcess, REALTIME_PRIORITY_CLASS);
start := now;
init;
makeprimes;
messageboxA(0,
PChar(Format('sum = %d, Time = %.2f',
[r, (now - start) * SecsPerDay])),
'prime', 0);
end.
yaos 2004-06-27
- 打赏
- 举报
也就是有4种情况
位数组，小cache
位数组，大cache
字节数组，小cache
字节数组，大cache
位数组，小cache
位数组，大cache
字节数组，小cache
字节数组，大cache
shines77 2004-06-25
- 打赏
- 举报
16进制？什么意思？不明白？
你先看下我的程序，我的循环块都是16进制的啊
最近2天生活有点乱，下了beta3，可服务器就上去了一下就当了，就看见了个开场动画
你先看下我的程序，我的循环块都是16进制的啊
最近2天生活有点乱，下了beta3，可服务器就上去了一下就当了，就看见了个开场动画
yaos 2004-06-25
- 打赏
- 举报
你有时间做做吧
我最近没时间调试程序了
另外，我们能不能从16进制考虑问题 ：）
我最近没时间调试程序了
另外，我们能不能从16进制考虑问题 ：）
shines77 2004-06-25
- 打赏
- 举报
今天查了一下，的确 BT 和 BTs (s为任一字母) 指令可以对位检测和赋值，应该有得优化
yaos 2004-06-24
- 打赏
- 举报
386汇编就可以根据偏移来存取位信息，才1个时钟，严重优化阿
yaos 2004-06-23
- 打赏
- 举报
嘿嘿，我还不打算上汇编，调试比较麻烦，暂时用Byte吧 ：）
yaos 2004-06-23
- 打赏
- 举报
多谢指教
确实是Bit压缩效率高，不过要牺牲点运算时间啊，除非用386 Bit运算指令
另外，昨天没想仔细，第一个块要做特殊处理，最好是直接数 ：）
素数表的生成也要特殊处理，打算利用筛的结果，需要一个初始素数表，大概
小于sqrt(CacheSize)就可以了，但是总的素数表至少要包含小于cacheSize的
所有素数！！
程序中，已经忽略了偶数了 ：）
j + = 2 * P[i]
至于你再次强调的Bit标志问题，还是用汇编做比较合适，性能提高很大的哦
确实是Bit压缩效率高，不过要牺牲点运算时间啊，除非用386 Bit运算指令
另外，昨天没想仔细，第一个块要做特殊处理，最好是直接数 ：）
素数表的生成也要特殊处理，打算利用筛的结果，需要一个初始素数表，大概
小于sqrt(CacheSize)就可以了，但是总的素数表至少要包含小于cacheSize的
所有素数！！
程序中，已经忽略了偶数了 ：）
j + = 2 * P[i]
至于你再次强调的Bit标志问题，还是用汇编做比较合适，性能提高很大的哦
shines77 2004-06-23
- 打赏
- 举报
WoW Beta2的私服我玩过了，打怪和不能打怪（看风景的），Beta3的私服也出来了，可惜硬盘最近有点紧张，暂时没有下载
其实你没有忽略偶数，你只是跳过，真正忽略应该都不保存它的标记位，这样你的CacheSize要大一倍，还有for中i++, i+=2的效率是不一样的
质数表是这样的，要求N内的质数，至少需要求出sqrt(N)内的所有质数，这个好办，对于2^32（即DWORD），只要求出2^16(65536)内的质数就可以了，所以我的循环块(3*5*7*4)×K=420×K（K为系数，使420*K刚好大于sqrt(N)，1<=K<=10），轻易就可以覆盖了这个65536内的质数，因为我的一个Byte覆盖了16个数的标记位（为什么看前面的bit压缩），所以420×K的循环块其实相当于420*16*K=6720*K，K=10时，＝67200，还比65536略多，所以除了前面一个420×K是特殊处理以外（这是一个小的循环节，重新看了下我的代码里是INIT_BLOCK_SIZE ((3*5*7 * 4) * 1)），真正开始循环我是从420*16×K+1开始的，循环块的大小是CACHE_BLOCK_SIZE (60060) = (3*5*7*11*4)*13，其实后来我改为73920稍快一些
至于Bit标志的问题，不用汇编也可以的，因为编译器基本上可以给你一个比较好（只能说还行，否则编译器也不是吃素的）的效率，当然我也还没有用汇编试试，多数的情况下没有很高的技巧或者明显汇编可以优化的地方，是不会有很大的效率优化
其实你没有忽略偶数，你只是跳过，真正忽略应该都不保存它的标记位，这样你的CacheSize要大一倍，还有for中i++, i+=2的效率是不一样的
质数表是这样的，要求N内的质数，至少需要求出sqrt(N)内的所有质数，这个好办，对于2^32（即DWORD），只要求出2^16(65536)内的质数就可以了，所以我的循环块(3*5*7*4)×K=420×K（K为系数，使420*K刚好大于sqrt(N)，1<=K<=10），轻易就可以覆盖了这个65536内的质数，因为我的一个Byte覆盖了16个数的标记位（为什么看前面的bit压缩），所以420×K的循环块其实相当于420*16*K=6720*K，K=10时，＝67200，还比65536略多，所以除了前面一个420×K是特殊处理以外（这是一个小的循环节，重新看了下我的代码里是INIT_BLOCK_SIZE ((3*5*7 * 4) * 1)），真正开始循环我是从420*16×K+1开始的，循环块的大小是CACHE_BLOCK_SIZE (60060) = (3*5*7*11*4)*13，其实后来我改为73920稍快一些
至于Bit标志的问题，不用汇编也可以的，因为编译器基本上可以给你一个比较好（只能说还行，否则编译器也不是吃素的）的效率，当然我也还没有用汇编试试，多数的情况下没有很高的技巧或者明显汇编可以优化的地方，是不会有很大的效率优化
shines77 2004-06-23
- 打赏
- 举报
更正，大体的提纲
shines77 2004-06-23
- 打赏
- 举报
如果你上面的代码只是大体的提供的话，还可以，不过你的Cache应该是这样定义（按照你的分块的逻辑）才对：
#define CacheSize (2 * 3 * 5 * 7 * 11 * 32)
还有，如果你要把第一个分块的位置放在从0开始的话，要注意第一个分块是特殊的，你注意到了吗？因为你把2,3,5,7,11都筛掉了（可是他们是质数），所以一般第一个分块不从0开始，或者做特别处理
#define CacheSize (2 * 3 * 5 * 7 * 11 * 32)
还有，如果你要把第一个分块的位置放在从0开始的话，要注意第一个分块是特殊的，你注意到了吗？因为你把2,3,5,7,11都筛掉了（可是他们是质数），所以一般第一个分块不从0开始，或者做特别处理
shines77 2004-06-22
- 打赏
- 举报
我只能说很遗憾，你是用一个char来表示一个标记位（就是一个BYTE），而为了压缩数据，可以使用bit来做标记位，这样不是省了7/8的存储空间了吗？再者，偶数中只有2是质数，我们完全可以一开始就把偶数给去掉，标记位中第一个bit存储的是1的标记位，依次是3,5,7,9,......，偶数根本没有必要处理，这样就又省了1/2，为什么要这么省呢？虽然这样读写标记位处理是复杂了，而节省的空间却可以使用较大的循环块，因而总体的速度还是快了很多
yaos 2004-06-22
- 打赏
- 举报
嘿嘿，刚放上去，就发现了一个错误 :)
n < 13的计算是错的，似乎也不能根据n < 13划分情况
n < 13的计算是错的，似乎也不能根据n < 13划分情况
yaos 2004-06-22
- 打赏
- 举报
#include <iostream.h>
#define CacheSize 4 * 16 * 3 * 5 * 7 * 11
unsigned char PCache[CacheSize];
unsigned char PBuffer[CacheSize];
unsigned long n;
unsigned long sum = 0;
unsigned long P[]; //素数表
//分块筛法计算n --> n + p - 1的素数
void BlockSieve(unsigned long n, unsigned long m)
{
unsigned long i, j, k, l, r, t;
for (i = 0; i < CacheSize; i ++)
PBuffer[i] = PCache[i];
k = sqrt(n);
l = sqrt(N + p - 1);
for (i = 5; P[i] <= k; i ++)
{
//标记P[i]的倍数, 在PBuffer中
r = (n / P[i]) * P[i];
if (r < n)
r = r + P[i];
if (r % 2 == 0)
r = r + P[i];
for (j = r - n; j <= m - 1; j + = 2 * P[i])
PBuffer[j] = 0;
}
for (; P[i] <= l; i ++)
{
//标记P[i]的倍数
for (j = P[i] * P[i] - n; j < m - 1; j += 2 * P[i])
PBuffer[j] = 0;
}
//输出
for (i = 0; i < CacheSize; i ++)
if (PBuffer[i])
{
//输出
sum ++;
}
}
void initPrime(unsigned long n)
{
unsigned long ps = n / ln(n) * 1.2; //大概素数个数
if (ps < 10) ps = 10;
P = new unsigned long[ps];
//同本程序的筛法
...
}
void initCache(void);
{
unsigned long i;
for (i = 0; i < CacheSize; i ++)
{
if ((i % 2) == 0)
PCache[i] = 0;
else
if ((i % 3) == 0)
PCache[i] = 0;
else
if ((i % 5) == 0)
PCache[i] = 0;
else
if ((i % 7) == 0)
PCache[i] = 0;
else
if ((i % 11) == 0)
PCache[i] = 0;
else
PCache[i] = 1;
}
}
void main(void)
{
unsigned i, j, m, k;
initCache();
cin >> n;
//初始化素数表
initPrime(sqrt(n));
if (n < 1)
cout << "0";
else
if (n < 13)
for (i = 0; i < n; i ++)
{
if (PCache[i])
{
//输出
sum ++;
}
}
else
{
m = n / (cacheSize);
k = n % (CacheSize);
if (m > 0)
for (i = 0; i < m; i ++)
BlockSieve(i * CacheSize, CacheSize);
BlockSieve(m * CacheSize, k);
}
//输出sum
cout >> endl >> sum;
}
先写提纲,后完善 :)
#define CacheSize 4 * 16 * 3 * 5 * 7 * 11
unsigned char PCache[CacheSize];
unsigned char PBuffer[CacheSize];
unsigned long n;
unsigned long sum = 0;
unsigned long P[]; //素数表
//分块筛法计算n --> n + p - 1的素数
void BlockSieve(unsigned long n, unsigned long m)
{
unsigned long i, j, k, l, r, t;
for (i = 0; i < CacheSize; i ++)
PBuffer[i] = PCache[i];
k = sqrt(n);
l = sqrt(N + p - 1);
for (i = 5; P[i] <= k; i ++)
{
//标记P[i]的倍数, 在PBuffer中
r = (n / P[i]) * P[i];
if (r < n)
r = r + P[i];
if (r % 2 == 0)
r = r + P[i];
for (j = r - n; j <= m - 1; j + = 2 * P[i])
PBuffer[j] = 0;
}
for (; P[i] <= l; i ++)
{
//标记P[i]的倍数
for (j = P[i] * P[i] - n; j < m - 1; j += 2 * P[i])
PBuffer[j] = 0;
}
//输出
for (i = 0; i < CacheSize; i ++)
if (PBuffer[i])
{
//输出
sum ++;
}
}
void initPrime(unsigned long n)
{
unsigned long ps = n / ln(n) * 1.2; //大概素数个数
if (ps < 10) ps = 10;
P = new unsigned long[ps];
//同本程序的筛法
...
}
void initCache(void);
{
unsigned long i;
for (i = 0; i < CacheSize; i ++)
{
if ((i % 2) == 0)
PCache[i] = 0;
else
if ((i % 3) == 0)
PCache[i] = 0;
else
if ((i % 5) == 0)
PCache[i] = 0;
else
if ((i % 7) == 0)
PCache[i] = 0;
else
if ((i % 11) == 0)
PCache[i] = 0;
else
PCache[i] = 1;
}
}
void main(void)
{
unsigned i, j, m, k;
initCache();
cin >> n;
//初始化素数表
initPrime(sqrt(n));
if (n < 1)
cout << "0";
else
if (n < 13)
for (i = 0; i < n; i ++)
{
if (PCache[i])
{
//输出
sum ++;
}
}
else
{
m = n / (cacheSize);
k = n % (CacheSize);
if (m > 0)
for (i = 0; i < m; i ++)
BlockSieve(i * CacheSize, CacheSize);
BlockSieve(m * CacheSize, k);
}
//输出sum
cout >> endl >> sum;
}
先写提纲,后完善 :)
加载更多回复（23）
