计算孪生素数对

tprime

孪生素数: 对于给定的素数p, 使p+2也为素数, 满足条件的 p, p+2 称为孪生素数(对) 详细介绍请查看 : http://baike.baidu.com/view/272607.htm 题目: 对于给定的k(k < 1000) 组整数对m, n (0 < m < n < 2^31) 求出每一组区间中的孪生素数对数, 即 PI2(m, n) 对于内存要求( < 100M), 时间要求( < 10s ) (版主的素数计算程序没有给出计算孪生素数, 希望提供此功能)

eMath

原帖由 tprime 于 2007-12-26 13:30 发表 (版主的素数计算程序没有给出计算孪生素数, 希望提供此功能)

我的素数计算程序 PrimeNumber.zip (18.65 KB, 下载次数: 24)

2007-12-27 07:56 上传

点击文件名下载附件

是我开发HugeCalc过程中的“副产品”， 由于发现网上有不少人苦于没有足够大的“素数表”而犯愁， 所以才决定把它独立出来与大家分享。 相对来说，孪生素数没有素数的适用面广，用户的需求不是很大（但部分用户可能非常迫切）， 所以就没有继续向这方面走，当然还有从程序的整体架构、输入输出界面等方面综合考量。 楼主提出“擂台”，最好是有个参比对象，比如一段代码、或一个程序等， 这样大家才能有目标，并在追逐极限过程中享受快乐，同时达到相互促进。

mathe

就我所知，直接产生孪生素数表并没有特别好的办法，可能其复杂度本来就同产生素数表相同。这可能也是为何通常我们看不到专门介绍产生孪生素数表的方法

tprime

//欢迎大家测试 // PI2[2^31] = 6810670 // AMD 3600+ 2.0G 1.6 s // intel PD 940 3.2G 1.8 s /***************************************************************** describition: input two number beg, end (0 < beg <= end < 2^32) claculate number of twin primes in the interval [beg, end] ******************************************************************/ # include # include # include # include # include # include # define MULTI_THREAD // multi-core CPU optimaztion # define THREAD_NUM 8 //the running threads # define NUM8 8 typedef unsigned int uint; typedef unsigned char uchar; typedef unsigned short ushort; # ifdef UCHAR typedef uchar utype; # define COMP 4 #define MASK 7 # else typedef ushort utype; # define COMP 5 #define MASK 15 # endif const int TBITSIZE = 1 << (COMP - 1); const int SQRTN = 6600; #define MASKN(n) (1 << ((n >> 1) & MASK)) //const uint DEFAULT_N = 10000000; const uint DEFAULT_N = (1u << 31) - 1; //the input range [0, DEFAULT_N] utype bits2[1 << TBITSIZE]; //num of twin primes int Prime[SQRTN]; //primes in interval[0, SQRTN] int sieve(int); int PI2(int, int); #ifdef S6 #define FACT 15015 * 14 // 3 * 5 * 7 * 11 * 13 = 15015 #define SP 6 #else #define FACT 255255 // 3 * 5 * 7 * 11 * 13 * 17 = 255255 #define SP 7 #endif #define BLOCKSIZE (FACT << 1) #define BASEMEM ((BLOCKSIZE >> COMP) + TBITSIZE) utype BaseTpl[BASEMEM]; //the table the fist SP primes is removed /*****************************************************************/ char tablep30[32]; char tabletp30[64]; char multp2[NUM8 * 6][6]; const uchar pri30[] = { 1, 7, 11, 13, 17, 19, 23, 29 }; const uchar twinp30[] = { 1, 11, 13, 17, 19, 29 }; //const char twinp60[] = {1, 11, 13, 17, 19, 29, 31, 41, 43, 47, 49}; #ifdef MULTI_THREAD struct ThreadInfo { int beg, end; int primenums; }Threadparam[THREAD_NUM * 2 + 2]; #ifdef _WIN32 # include DWORD WINAPI Win32ThreadFunc(LPVOID pinfo) #else # include void* POSIXThreadFunc(void *pinfo) #endif { ThreadInfo *pThreadInfo = (ThreadInfo *)(pinfo); pThreadInfo->primenums = PI2(pThreadInfo->beg, pThreadInfo->end); // printf("Thread: PI2[%10d, %10d] = %d\n", pThreadInfo->beg, pThreadInfo->end, pThreadInfo->primenums); return 0; } int multiThread(int threadnums, int beg, uint end) { int i, primenums = 0; Threadparam[0].beg = beg; int tsize = (end - beg) / threadnums; for (i = 1; i < threadnums; i++) Threadparam[i].beg = Threadparam[i - 1].end = Threadparam[i - 1].beg + tsize; Threadparam[threadnums - 1].end = end; #ifdef _WIN32 HANDLE thredHand[THREAD_NUM * 2]; DWORD threadID[THREAD_NUM * 2]; for (i = 0; i < threadnums; i++){ thredHand[i] = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE)Win32ThreadFunc, (LPVOID)(&Threadparam[i]), 0, &threadID[i]); if (thredHand[i] == NULL) printf("create Win32 thread error = %d\n", GetLastError()); } WaitForMultipleObjects(threadnums, thredHand, true, INFINITE); for (i = 0; i < threadnums; i++){ primenums += Threadparam[i].primenums; CloseHandle(thredHand[i]); } #else pthread_t pid[THREAD_NUM * 2]; for (i = 0; i < threadnums; i++){ int error = pthread_create(&pid[i], NULL, POSIXThreadFunc, &Threadparam[i]); if ( error != 0 ) printf("create pthread error = %d\n", error); } for (i = 0; i < threadnums; i++){ pthread_join(pid[i], NULL); primenums += Threadparam[i].primenums; } #endif return primenums; } #endif void initmulp2(int t30i, int p30i) { int pgap = 0, lastgap = 0; for (int i = 0, mod30 = 0; i < sizeof(twinp30) / sizeof(twinp30[0]); ){ if (tabletp30[mod30] > 0 && pgap != lastgap){ multp2[t30i * NUM8 + p30i][i++] = pgap - lastgap; lastgap = pgap; } pgap += 2; mod30 = (twinp30[t30i] + pgap * pri30[p30i]) % 30; } } void inline initMask(utype mask[], int start, uint len) { int srid = start % BLOCKSIZE; memcpy(mask, BaseTpl, sizeof(BaseTpl)); if (start < 32) *((ushort *)&mask[0]) = ~0xcb6e; if (srid) mask[srid >> COMP] |= (MASKN(srid) - 1); if (len != BLOCKSIZE) mask[len >> COMP] |= ~(MASKN(len) - 1); } int piRange2(int start, uint len = BLOCKSIZE) { static char prime2[BASEMEM]; int srid = start % BLOCKSIZE; len += srid; uint next, primenums = 0; const bool ok = start >= BLOCKSIZE; const int sqrtn = (int)sqrt((float)(start + len)) + 1; utype mask[BASEMEM + 1]; initMask(mask, start--, len); for (int i = SP + 1, beg = srid - 1, p = Prime[i]; p < sqrtn; p = Prime[++i]){ if (ok){ next = beg + p - start % p; if ((next & 1) == 0) next += p; } else next = p * p; int buf[8]; const uchar nexti[] = { 1, 2, 3, 4, 5, 0}; while ( tabletp30[next % 30] == 0 ) next += p * 2; int st = 5, ps = tablep30[p % 30] - 1 + (tabletp30[next % 30] - 1) * 8; for (int k = 0; k < 6; k++) buf[k] = p * multp2[ps][k]; for (; next < len; next += buf[st = nexti[st]]) mask[next >> COMP] |= MASKN(next); } int size = len >> COMP; for (int k = (srid >> COMP); k <= size; k++){ primenums += bits2[mask[k]]; if (mask[k] < (1 << (TBITSIZE - 1)) && !(mask[k + 1] & 1)) primenums++; } int blocki = ++start / BLOCKSIZE; if (!srid && !(mask[0] & 1)){ prime2[blocki] |= 1; if (prime2[blocki - 1] > 1 /*&& blocki > 0*/) primenums++; } if (len == BLOCKSIZE && !(mask[(BLOCKSIZE >> COMP)] & MASKN(BLOCKSIZE - 1))) prime2[blocki] |= 2; return primenums; } int PI2(int beg, int end) { int primenums = 0; if (beg > end || end > DEFAULT_N) return -1; if (beg / BLOCKSIZE == end / BLOCKSIZE){ primenums += piRange2(beg, end - beg + 1); return primenums; } int size = end % BLOCKSIZE; if (size != 0) primenums += piRange2(end - size, size + 1); size = beg % BLOCKSIZE; if (size != 0){ primenums += piRange2(beg, BLOCKSIZE - size); beg += BLOCKSIZE - size; } beg /= BLOCKSIZE, end /= BLOCKSIZE; for (int i = beg; i < end; i++){ static int pcache[DEFAULT_N / BLOCKSIZE + 2]; if (pcache[i] == 0) pcache[i] = piRange2(i * BLOCKSIZE); primenums += pcache[i]; } return primenums; } void InitTable( ) { int i = 0; int bitsize = (1 << TBITSIZE); //init bits2 for (i = 1; i < bitsize; i++){ uint mask2 = 3; while (mask2 < (1 << TBITSIZE)){ if ((i & mask2) == 0) bits2[i]++; mask2 <<= 1; } } //init tablep30 for (i = 0; i < 8; i++) tablep30[pri30[i] % 30] = i + 1; //init twinp30 for (i = 0; i < 6; i++) tabletp30[twinp30[i] % 30] = i + 1; //init multp2 for (int m = 0; m < 6; m++) for (int r = 0; r < 8; r++) initmulp2(m, r); } int test(int tp, int t = 1000) { clock_t tstart = clock(); srand(time(0)); for (int i = 1; i <= t; i++){ int beg = rand() * rand() % (DEFAULT_N); int end = rand() * rand() % (DEFAULT_N); if (beg > end){ int t = beg; beg = end; end = t; } int p2 = PI2(beg, end); printf("case%4d : PI2(%9d %9d) = %d\n",i, beg, end, p2); } printf("test case number %d, total time use %d ms\n", t, clock() - tstart); //for test input data #ifdef FILEINPUT freopen("testr.txt", "r", stdin); freopen("testw.txt", "w", stdout); int beg, end; while (scanf ("%d %d", &beg, &end) == 2 && beg <= end){ tstart = clock(); int primenums = PI2(beg, end); printf("PI2[%9d, %9d] = ", beg, end); printf("%9d, time use %ld ms\n", primenums, clock() - tstart); } #endif return 0; } int main(int arg, char *argc[]) { clock_t tstart = clock(); int primenums = 0; int threadnum = THREAD_NUM; InitTable(); sieve(DEFAULT_N); #ifdef MULTI_THREAD if (arg > 1) threadnum = atoi(argc[1]); if (DEFAULT_N / (2 * BLOCKSIZE) > threadnum) primenums = multiThread(threadnum, 0, DEFAULT_N); else primenums = PI2(0, DEFAULT_N); #else primenums = PI2(0, DEFAULT_N); #endif printf("PI2[0 - %u] : primes = %d, init table time use %ld ms\n", (uint)DEFAULT_N, primenums, clock() - tstart); test(threadnum); //for test the result return 0; } int sieve(int maxn) { uint p, primes = 1; uint maxp = (uint)sqrt((float)maxn) + 19; Prime[1] = 2; for (p = 3; p < maxp; p += 2){ if ( !(BaseTpl[p >> COMP] & MASKN(p)) ){ Prime[++primes] = p; for (uint j = p * p; j < maxp; j += p << 1) BaseTpl[j >> COMP] |= MASKN(j); } } memset(BaseTpl, 0, sizeof(BaseTpl)); for (int i = 2; i <= SP; i++){ for (int j = Prime[i], p = j; j < BLOCKSIZE; j += p << 1) BaseTpl[j >> COMP] |= MASKN(j); } Prime[primes + 1] = maxn; Prime[0] = primes; BaseTpl[BLOCKSIZE >> COMP] |= ~(MASKN(BLOCKSIZE) - 1); BaseTpl[(BLOCKSIZE >> COMP) + 1] = ~0; return primes; }

admin

刚才发现你的帖子后半部是斜体，估计是“[i]”被误认为成Discuz!代码了， 我给你编辑了下，只要发帖时在“禁用 Discuz!代码”前打钩即可。

mathe

原来你是只需要计算数目，而不需要找出所有的素数对

tprime

找出来实现也很容易, 位数组已经筛掉不满足条件的，考虑到需要较多内存就没实现.

tprime

原帖有一些小bug, 尤其是多线程时会少算1, 进几天改进了算法, 如果只是计数PI2(n,) 速度提高100% 测试结果: Intel PD 940 PI2(2^31) 650ms, AMD 3600+ PI2(2^31) 500ms 下一个版本将实现存储, 有可能的话再加UI。 计算PI(x)有很多较快的算法, 写的好的程序1秒可以算出 PI(10^12),但计算PI2(x)不知道是否有比筛选法更好的？

gxqcn

我见过的计算pi(x)速度最快的是这篇论文： computing-pi-x-the.pdf (189.21 KB, 下载次数: 28)

2008-1-3 10:27 上传

点击文件名下载附件
快速计算pi(x)算法

COMPUTING pi((x)): THE MEISSEL-LEHMER METHOD 作者为：J. C. Lagarias(AT&T Bell Laboratories,Murray Hill, New Jersey); 　　　　V. S. Miller(IBM Watson Research Center,Yorktown Heights, New York); 　　　　A. M. Odlyzko(AT&T Bell Laboratories,Murray Hill, New Jersey) 关于PI2(x)，也许正如 mathe 所说的“并没有特别好的办法”，至少我没有看到太高效的算法。 tprime，如果你愿意，在你的程序定型后，可以考虑加入到数学研发网原创软件中进行统一宣传，因为我的收费主页空间还有些剩余（注意：因空间限制，只会收录一些优秀的、精巧的原创作品）。

tprime

再次更新的算法, 在Intel PD 940 3.2G 计算PI2(10^12), 需要6 分钟, 一个晚上估计能算到PI2(10^14) 目前国外采用分布式算法能算到PI2(10^17), HaHa! 还有很多工作要做: 计算n生素数比(n-1)生素数快2-3倍, (n > 2 && n < 7)