极大素数间隔问题

liangbch · 发表于 2008-3-19 15:31:34

不妨定义这样的需求，或许更有用些。对于每一个偶数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

gxqcn · 发表于 2008-3-19 19:32:09

原帖由 liangbch 于 2008-3-19 15:31 发表不妨定义这样的需求，或许更有用些。对于每一个偶数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 尽可能小”更恰当，与其它参考资料比较相融。

无心人 · 发表于 2008-3-20 09:10:25

一个模子扣出来的题目并没有区别啊也不见得更容易假设求k=2000的，如何操作？肯定在50位以内， 32位以上

liangbch · 发表于 2008-3-20 11:41:27

原帖由 无心人 于 2008-3-20 09:10 发表一个模子扣出来的题目并没有区别啊也不见得更容易

我没有仔细看题目，实际上我的需求和你的是一致的。

无心人 · 发表于 2008-3-21 09:24:30

program p27;
//发现在10^14以内的极大素数间隔
type
primearray = array[1..664800] of longword;
primesieve = array[1..4 * 3 * 5 * 7 * 11 * 13 * 17 * 19] 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[1..8] of integer;
begin
pp[1] := 2;
pp[2] := 3;
pp[3] := 5;
pp[4] := 7;
pp[5] := 11;
pp[6] := 13;
pp[7] := 17;
pp[8] := 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[j] = 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 [eax], ebx
mov dword ptr [eax + 4], edx
pop ebx
end;
end;
procedure findsmallprime;
var
i, j, k: longword;
begin
sp[1] := 3;
sp[2] := 5;
sp[3] := 7;
sp[4] := 11;
csp := 4;
k := 13;
while (k <= 9999999) do
begin
j := 1;
i := round(sqrt(k) + 0.5);
while sp[j] <= i do
begin
if k mod sp[j] = 0 then break;
inc(j);
end;
if k mod sp[j] <> 0 then
begin
inc(csp);
sp[csp] := 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 [stop]
fsqrt
fadd qword ptr [c01]
fistp dword ptr [m]
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 [j + 4]
mov ecx, tmp2
div ecx
mov ebx, eax
mov eax, dword ptr [j]
div ecx
mul ecx
mov dword ptr [k], eax
mov eax, ebx
mov ebx, edx
mul ecx
add eax, ebx
mov dword ptr [k + 4], 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[k- start + 1] <> 0 then
p[k - start + 1] := 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.

复制代码

mathe · 发表于 2008-3-21 09:41:16

17071454401 = 880*19399380+1 为什么要从这里开始？有什么先验知识在这里？还有为什么选择4 * 3 * 5 * 7 * 11 * 13 * 17 * 19，其实应该 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19就可以了呀？

无心人 · 发表于 2008-3-21 09:58:20

哈因为是之前的已计算了啊至于选4 * 3 * 5 * 7 * 11 * 13 * 17 * 19因为是32位整数的倍数是02年4月的东西，现在肯定不这么写

mathe · 发表于 2008-3-24 11:03:34

原帖由 无心人 于 2008-3-19 13:55 发表我能调度但做仔细的编码 1、没相关算法 2、手里的汇编优化资料全E文的，理解有困难 3、通过这几个星期的练习汇编功力还有待加强啊

这个算法不是问题，袁本身的算法他大概描述过，其实不是很复杂。我们可以认为分成两个部分，第一部分是搜索过程，对于每个给定的低N/4位，可以唯一确定开平方后的结果（基本唯一确定）。对于这部分代码，袁的算法不算最优，我们可以继续优化到基本上不出现乘法运算。但是，在找到开平方后的结果，我们还需要验证这个结果是否正确的，而这部分验证过程的时间复杂度和上面过程是相同的，但是其中乘法计算（相当于一个比较大的整数的平方计算，基本上是20到30位的10进制数的平方，最好有快速的10进制计算）我没法优化掉，所以这一部分的时间很难优化掉。也正因为这样，我觉得要远远超越他的结果很难，但是考虑到你们可以对乘法做很好的优化的话（不知道比普通C代码写出来可以快多少倍），那么应该可以平均速度比他快一些。但是在算法复杂度上，并没有明显的改变。

无心人 · 发表于 2008-3-24 16:02:16

如果测算乘法 8位一组极限情况下，就是32位4组乘了如果有好算法能快速从 x * (2^64 / 10^8)中得到 x / 10^8 商，那10进乘法并不困难，否则可能时间是二进乘法的2倍甚至更多 ======================================= 或者做6位一组的乘，这个纠正系数是10995很小了

mathe · 发表于 2008-3-24 16:21:00

这个div by constant转化为乘法总是可以做的，你可以让C编译器编译一下，看优化以后的汇编代码就应该能够得到

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 极大素数间隔问题

可爱的程序终于在CSDN搜索到了：）

[讨论] 极大素数间隔问题

可爱的程序终于在CSDN搜索到了 ：）

可爱的程序终于在CSDN搜索到了：）