对于每一个偶数K, 求出相邻的质数对(p1,p2),使得p2-p1=K 且 p2 尽可能小。
如对于K=2,p1=3, p2=5
K=4,p1=7, p2=11
K=14, p1=113, p2=127 原帖由 liangbch 于 2008-3-19 15:31 发表 http://images.5d6d.net/dz60/common/back.gif
不妨定义这样的需求,或许更有用些。
对于每一个偶数K, 求出相邻的质数对(p1,p2),使得p2-p1=K 且 p2 尽可能小。
如对于K=2,p1=3, p2=5
K=4,p1=7, p2=11
K=14, p1=113, p2=127
感觉改成使“p1 尽可能小”更恰当,与其它参考资料比较相融。 一个模子扣出来的题目
并没有区别啊
也不见得更容易
假设求k=2000的,如何操作?
肯定在50位以内, 32位以上 原帖由 无心人 于 2008-3-20 09:10 发表 http://images.5d6d.net/dz60/common/back.gif
一个模子扣出来的题目
并没有区别啊
也不见得更容易
我没有仔细看题目,实际上我的需求和你的是一致的。
可爱的程序终于在CSDN搜索到了 :)
program p27;//发现在10^14以内的极大素数间隔
type
primearray = array of longword;
primesieve = array of byte;
const
c01: double = 0.5;
var
p: primesieve;
t: primesieve;
sp: primearray;
lp, rp: int64;
start, stop: int64;
prev: int64;
pf: TextFile;
s, csp: integer;
cpuA, cpuB: int64;
procedure compsieve;
var
i, j: integer;
pp: array of integer;
begin
pp := 2;
pp := 3;
pp := 5;
pp := 7;
pp := 11;
pp := 13;
pp := 17;
pp := 19;
for i := 1 to 4 * 3 * 5 * 7 * 11 * 13 * 17 * 19 do
begin
t[ i ] := 1;
for j := 1 to 8 do
if i mod pp = 0 then
begin
t[ i ] := 0;
end;
end;
end;
procedure comptemp;
var
tf: file of primesieve;
begin
assignfile(tf, 'psieve.dat');
{\$I-}
reset(tf);
read(tf, t);
close(tf);
{\$I+}
if IOResult <> 0 then
begin
compsieve;
rewrite(tf);
write(tf, t);
closefile(tf);
end;
end;
procedure loadtemp;
begin
move(t, p, sizeof(p));
end;
function IsPrime(TestPrime: int64): boolean;
var
i: integer;
k: int64;
begin
i := 1;
k := round(sqrt(TestPrime * 1.0) + 0.5);
while sp[ i ] <= k do
begin
if TestPrime mod sp[ i ] = 0 then
begin
IsPrime := False;
exit;
end;
inc(i);
end;
IsPrime := True;
end;
procedure getcputime(var cputime: int64);
begin
asm
push ebx
push eax
db \$0f, \$31
mov ebx, eax
pop eax
mov dword ptr , ebx
mov dword ptr , edx
pop ebx
end;
end;
procedure findsmallprime;
var
i, j, k: longword;
begin
sp := 3;
sp := 5;
sp := 7;
sp := 11;
csp := 4;
k := 13;
while (k <= 9999999) do
begin
j := 1;
i := round(sqrt(k) + 0.5);
while sp <= i do
begin
if k mod sp = 0 then break;
inc(j);
end;
if k mod sp <> 0 then
begin
inc(csp);
sp := k;
end;
k := k + 2;
end;
end;
procedure loadprime;
var
p: file of primearray;
begin
assignfile(p, 'prime.dat');
{\$I-}
reset(p);
read(p, sp);
{\$I+}
if IOResult <> 0 then
begin
FindSmallPrime;
rewrite(p);
write(p, sp);
end;
csp := 664578;
closefile(p);
end;
procedure sieve;
var
j, k, n: int64;
i, l: integer;
m: longword;
tmp1: longword;
tmp2: longword;
begin
stop := start + 4 * 3 * 5 * 7 * 11 * 13 * 17 * 19 - 1;
// getcputime(cpuA);
loadtemp;
//m := round(sqrt(stop * 1.0)+0.5)
asm
fild qword ptr
fsqrt
fadd qword ptr
fistp dword ptr
end;
i := 8;
repeat
tmp2 := sp[ i ];
tmp1 := tmp2 * 2;
n := tmp2 * tmp2;
if n < start then j := start else j := n;
// k := (j div tmp2) * tmp2;
asm
push ebx
push ecx
mov edx, 0
mov eax, dword ptr
mov ecx, tmp2
div ecx
mov ebx, eax
mov eax, dword ptr
div ecx
mul ecx
mov dword ptr , eax
mov eax, ebx
mov ebx, edx
mul ecx
add eax, ebx
mov dword ptr , eax
pop ecx
pop ebx
end;
if k < start then k := k + tmp2;
if k mod 2 = 0 then k := k + tmp2;
while (k <= stop) do
begin
if p <> 0 then
p := 0;
k := k + tmp1;
end;
inc(i);
until (sp[ i ] > m) or (i > csp);
// getcputime(cpuB);
// writeln(' Sieve cpu code time:', cpuB - cpuA);
l := prev - start + 1;
for i := 1 to 4 * 3 * 5 * 7 * 11 * 13 * 17 * 19 do
begin
if p[ i ] = 1 then
begin
if i - l >= s then
begin
lp := start + l - 1; rp := start + i - 1; s := i - l;
writeln;
writeln(lp:14, rp:14, s:14);
writeln(pf, lp:14, rp:14, s:14);
flush(pf);
end;
l := i;
end;
end;
prev := start + l - 1;
start := stop + 1;
end;
var
bf: TextFile;
bk: longword;
begin
loadprime;
// writeln('加载素数表成功!');
comptemp;
// writeln('准备模版成功!');
bk := 1;
assignfile(bf, 'prime.ini');
{\$I-}
reset(bf);
readln(bf, Start);
if start mod 19399380 <> 1 then
start := (start div 19399380) * 19399380 + 1;
readln(bf, s);
closefile(bf);
{\$I+}
if IOResult <> 0 then
begin
start := int64(17071454401);
s := 382;
end;
Prev := Start;
while not IsPrime(Prev) do
Prev := Prev - 2;
// writeln('加载开始点成功,start = ', start:14, ', prev = ', prev:14);
assignfile(pf, 'p2.txt');
append(pf);
repeat
sieve;
// write('.');
// writeln('现在进行下一个循环! prev = ', prev);
inc(bk);
if bk >= 100 then
begin
bk := 1;
rewrite(bf);
writeln(bf, start);
writeln(bf, s);
closefile(bf);
end;
until prev >= int64(100000000000000);
closefile(pf);
end. 17071454401 = 880*19399380+1
为什么要从这里开始?有什么先验知识在这里?
还有为什么选择4 * 3 * 5 * 7 * 11 * 13 * 17 * 19,其实应该
2 * 3 * 5 * 7 * 11 * 13 * 17 * 19就可以了呀? 哈 因为是之前的已计算了啊
至于选4 * 3 * 5 * 7 * 11 * 13 * 17 * 19因为是32位整数的倍数
是02年4月的东西, 现在肯定不这么写 原帖由 无心人 于 2008-3-19 13:55 发表 http://images.5d6d.net/dz60/common/back.gif
我能调度
但做仔细的编码
1、没相关算法
2、手里的汇编优化资料全E文的, 理解有困难
3、通过这几个星期的练习汇编功力还有待加强啊
这个算法不是问题,袁本身的算法他大概描述过,其实不是很复杂。
我们可以认为分成两个部分,第一部分是搜索过程,对于每个给定的低N/4位,可以唯一确定开平方后的结果(基本唯一确定)。对于这部分代码,袁的算法不算最优,我们可以继续优化到基本上不出现乘法运算。
但是,在找到开平方后的结果,我们还需要验证这个结果是否正确的,而这部分验证过程的时间复杂度和上面过程是相同的,但是其中乘法计算(相当于一个比较大的整数的平方计算,基本上是20到30位的10进制数的平方,最好有快速的10进制计算)我没法优化掉,所以这一部分的时间很难优化掉。也正因为这样,我觉得要远远超越他的结果很难,但是考虑到你们可以对乘法做很好的优化的话(不知道比普通C代码写出来可以快多少倍),那么应该可以平均速度比他快一些。
但是在算法复杂度上,并没有明显的改变。 如果测算乘法
8位一组
极限情况下,就是32位4组乘了
如果有好算法能快速从 x * (2^64 / 10^8)中得到 x / 10^8 商,那10进乘法并不困难,否则可能时间是二进乘法的2倍甚至更多
=======================================
或者做6位一组的乘,这个纠正系数是10995很小了 这个div by constant转化为乘法总是可以做的,你可以让C编译器编译一下,看优化以后的汇编代码就应该能够得到