孪生素数的计算

无心人 · 发表于 2008-7-17 08:40:36

谢谢我在写一个东西准备挂服务器上算2000－3000位的四生素数呵呵

无心人 · 发表于 2008-7-17 10:24:40

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "gmp.h"
unsigned long primes[1280]; //小于10000素数，实际上有1229个
unsigned long quadprimes[4800]; //小于10^8四生素数，实际上有4768组
unsigned long sieve[30030];
unsigned long template[30030];
void
InitPrimes (void)
{
primes[0] = 2;
primes[1] = 3;
primes[2] = 5;
primes[3] = 7;
primes[4] = 11;
primes[5] = 13;
primes[6] = 17;
primes[7] = 19;
primes[8] = 23;
primes[9] = 29;
primes[10] = 31;
primes[11] = 37;
primes[12] = 41;
primes[13] = 43;
primes[14] = 47;
primes[15] = 53;
primes[16] = 59;
primes[17] = 61;
primes[18] = 67;
primes[19] = 71;
primes[20] = 73;
primes[21] = 79;
primes[22] = 83;
primes[23] = 89;
primes[24] = 97;
}
void
InitSieve (void)
{
unsigned long i;
for (i = 1; i <= 9999; i += 2)
sieve[i] = 1;
}
void
Sieve (void)
{
unsigned long i, j, start;
for (i = 0; i <= 24; i++)
{
start = primes[i] * primes[i];
for (j = start; j <= 9999; j += 2 * primes[i])
sieve[j] = 0;
}
}
void
StoreSieveResult (void)
{
unsigned long i, k;
k = 25;
for (i = 101; i <= 9999; i += 2)
if (sieve[i])
primes[k++] = i;
for (; k < 1280; k++)
primes[k] = 0;
primes[1229] = 10007; //以这个素数结束，否则下面的四生素数筛在接近10^8会超越1229素数的范围
}
void
Init (void)
{
InitPrimes ();
InitSieve ();
Sieve ();
StoreSieveResult ();
}
void
InitQuadPrimes (void)
{
unsigned long i, k = 0;
for (i = 0; i < 1229; i++)
if ((primes[i + 1] - primes[i] == 2) && (primes[i + 2] - primes[i] == 6)
&& (primes[i + 3] - primes[i] == 8))
quadprimes[k++] = primes[i]; //共12组
}
void
InitTemplate (void)
{
unsigned long i;
for (i = 1; i < 30030; i += 2)
template[i] = 1;
for (i = 0; i < 30030; i += 2)
template[i] = 0;
for (i = 3; i < 30030; i += 6)
template[i] = 0;
for (i = 5; i < 30030; i += 10)
template[i] = 0;
for (i = 7; i < 30030; i += 14)
template[i] = 0;
for (i = 11; i < 30030; i += 22)
template[i] = 0;
for (i = 13; i < 30030; i += 26)
template[i] = 0;
}
void
SearchQuadSieve (void)
{
unsigned long start, end;
unsigned long i, j, k, m, s_start;
InitTemplate ();
//首次筛过程
start = 0;
end = 30030;
memcpy (sieve, template, 30030 * sizeof (unsigned long));
for (j = 6; primes[j] * primes[j] <= end; j++)
{
s_start = primes[j] * primes[j];
for (k = s_start; k < 30030; k += 2 * primes[j])
sieve[k] = 0;
}
m = 12;
for (i = 10001; i < 30030; i += 10)
if ((sieve[i]) && (sieve[i + 2]) && (sieve[i + 6]) && (sieve[i + 8]))
quadprimes[m++] = i;
//完整筛过程，经测试剩余的数字是99999900-99999999，不存在四生素数
for (i = 1; i < 100000000 / 30030; i++)
{
start = i * 30030;
end = start + 30030;
memcpy (sieve, template, 30030 * sizeof (unsigned long));
for (j = 6; primes[j] * primes[j] <= start; j++)
{
s_start = (start / primes[j]) * primes[j] + primes[j];
if (s_start % 2 == 0)
s_start += primes[j];
s_start -= start;
for (k = s_start; k < 30030; k += 2 * primes[j])
sieve[k] = 0;
}
for (; primes[j] * primes[j] < end; j++)
{
s_start = primes[j] * primes[j] - start;
for (k = s_start; k < 30030; k += 2 * primes[j])
sieve[k] = 0;
}
s_start = (start / 30) * 30 + 11;
if (s_start < start)
s_start += 30;
s_start -= start;
for (k = s_start; k < 30030; k += 30)
if ((sieve[k]) && (sieve[k + 2]) && (sieve[k + 6]) && (sieve[k + 8]))
quadprimes[m++] = start + k;
}
}
void
product (mpz_t r, unsigned long p[], unsigned long max)
{
unsigned long i;
mpz_set_ui (r, 1);
if (max < 2)
return;
for (i = 0; primes[i] <= max; i++)
mpz_mul_ui (r, r, p[i]);
}
int
main (void)
{
unsigned long a, c;
mpz_t b, p;
Init ();
InitQuadPrimes ();
SearchQuadSieve ();
for (a = 0; a < 4800; a++)
if (quadprimes[a])
printf ("%6u ", quadprimes[a]);
printf ("\n");
return 0;
}

复制代码

无心人 · 发表于 2008-7-17 10:27:08

上面是俺写大四生素数到初始化部分的程序初始化了一个10000内素数表，一个100000000内四生素数表下面是四生素数的结果，不知道对不对

quadprimes.txt.bz2 (16.21 KB, 下载次数: 4)

无心人 · 发表于 2008-7-17 11:09:04

不包括输出是570ms P4 2.0G/1280M/ubuntu 8.04/GCC 4.2.3

tprime · 发表于 2008-7-17 12:30:58

这里有10亿的

tprime · 发表于 2008-7-17 12:32:26

眼馋你的服务器, 多少个节点? CPU配置?

无心人 · 发表于 2008-7-17 15:00:54

一个破服务器啊 XEON 2.0 *2/1024M 就是可24小时开机而已谢谢你的文档

tprime · 发表于 2008-12-23 16:35:13

终于完成计算[0-2^32]任意两数之间的孪生素数(比最快的版本慢1倍, 不同算法) 随机测试100万组数据还没发现错误.目前还不能支持孪生素数列表(输出太慢不打算做了) 欢迎大家下载玩并提出意见, 下面是截屏输出和自己的测试数据. sieve cache size = 2 * 30 * 7 * 11 * 13 * 17 * 12 / 60 = 199K PI2(2147483647) = 6810670, init twin prime table time use 1767.14 ms please input two number: (0<= beg <= end < 2^32): 1 2147483647 case 1 : PI2[1,2147483647] = 6810670, time use 3 ms 1234567890 2147483647 case 2 : PI2[1234567890,2147483647] = 2672471, time use 6 ms 123456789 987654321 case 3 : PI2[123456789,987654321] = 2856051, time use 5 ms

独木星空 · 发表于 2010-8-29 23:04:36

8楼的数据很有用

独木星空 · 发表于 2010-8-29 23:08:12

N值→→→→4(8)素数个数→孪生素数对数目→2.38128115124→误差→误差率CP 10000000→898→→→→→58755→→→→822→→→→→76→0.084632517t 20000000→1467→→→→→107246→→→→1369→→→→→98→0.0668029991 30000000→1951→→→→→152790→→→→1853→→→→→98→0.050230651cBW:& 40000000→2403→→→→→196566→→→→2300→→→→→103→0.042863088l;\`5/q 50000000→2846→→→→→239094→→→→2722→→→→→124→0.043569923{ 60000000→3257→→→→→280666→→→→3126→→→→→131→0.0402210625h 70000000→3646→→→→→321468→→→→3515→→→→→131→0.035929786~# 80000000→4033→→→→→361627→→→→3892→→→→→141→0.034961567Hm9CBn 90000000→4401→→→→→401236→→→→4259→→→→→142→0.0322653948V_I 100000000→4767→→→→→440366→→→→4617→→→→→150→0.031466331nL[QTQ 110000000→5115→→→→→479074→→→→4968→→→→→147→0.028739003X< 120000000→5441→→→→→517402→→→→5312→→→→→129→0.0237088770*L 130000000→5797→→→→→555388→→→→5650→→→→→147→0.025357944Wz( 140000000→6112→→→→→593062→→→→5982→→→→→130→0.021269634y\$N0Xt 150000000→6450→→→→→630450→→→→6309→→→→→141→0.021860465Gvh,uS 160000000→6792→→→→→667574→→→→6632→→→→→160→0.023557126~0l 170000000→7113→→→→→704453→→→→6951→→→→→162→0.0227752]OX' 180000000→7446→→→→→741104→→→→7266→→→→→180→0.0241740531]4 190000000→7794→→→→→777541→→→→7577→→→→→217→0.027841936Sgq 200000000→8096→→→→→813777→→→→7884→→→→→212→0.026185771Tv2n 210000000→8400→→→→→849824→→→→8189→→→→→211→0.025119048HEt 220000000→8699→→→→→885692→→→→8490→→→→→209→0.02402575Md&q 230000000→8978→→→→→921391→→→→8789→→→→→189→0.021051459 M 240000000→9270→→→→→956928→→→→9085→→→→→185→0.01995685 h3Dh 250000000→9565→→→→→992313→→→→9379→→→→→186→0.019445896dyYRqA 260000000→9836→→→→→1027551→→→→9670→→→→→166→0.016876779FOTa 270000000→10135→→→→→1062650→→→→9959→→→→→176→0.017365565T 280000000→10431→→→→→1097615→→→→10245→→→→→186→0.017831464W,"Mt 290000000→10701→→→→→1132452→→→→10530→→→→→171→0.015979815} 300000000→10972→→→→→1167166→→→→10813→→→→→159→0.014491433Q 310000000→11280→→→→→1201762→→→→11093→→→→→187→0.016578014=j 320000000→11589→→→→→1236243→→→→11372→→→→→217→0.018724653:gJ1 330000000→11862→→→→→1270615→→→→11649→→→→→213→0.0179565C6zb 340000000→12126→→→→→1304880→→→→11925→→→→→201→0.016575952;#nw+$ 350000000→12370→→→→→1339043→→→→12199→→→→→171→0.013823767lS 360000000→12632→→→→→1373107→→→→12471→→→→→161→0.012745408yg%s} 370000000→12900→→→→→1407075→→→→12742→→→→→158→0.0122480628 380000000→13164→→→→→1440950→→→→13011→→→→→153→0.011622607Zj|^Cw 390000000→13438→→→→→1474735→→→→13279→→→→→159→0.011832118l#]9A 400000000→13712→→→→→1508433→→→→13545→→→→→167→0.012179113[a3Lz 410000000→13957→→→→→1542045→→→→13810→→→→→147→0.010532349q 420000000→14247→→→→→1575575→→→→14074→→→→→173→0.012142907VWu0; 430000000→14516→→→→→1609025→→→→14337→→→→→179→0.012331221Z^|(x 440000000→14770→→→→→1642396→→→→14598→→→→→172→0.011645227j 450000000→15030→→→→→1675691→→→→14858→→→→→172→0.011443779| 460000000→15289→→→→→1708912→→→→15117→→→→→172→0.0112499181 470000000→15559→→→→→1742061→→→→15375→→→→→184→0.011825953b3:<2& 480000000→15823→→→→→1775139→→→→15632→→→→→191→0.012071036B|wI,j 490000000→16093→→→→→1808148→→→→15888→→→→→205→0.012738458u,?Pb2 500000000→16330→→→→→1841090→→→→16143→→→→→187→0.011451317A1 510000000→16579→→→→→1873967→→→→16396→→→→→183→0.01103806Sk6!kK 520000000→16816→→→→→1906779→→→→16649→→→→→167→0.009931018u}r\$b 530000000→17069→→→→→1939528→→→→16901→→→→→168→0.009842404m8Z8?

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

[讨论] 孪生素数的计算