计算孪生素数对
孪生素数: 对于给定的素数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 )
(版主的素数计算程序没有给出计算孪生素数, 希望提供此功能) 原帖由 tprime 于 2007-12-26 13:30 发表 http://emath.5d6d.com/images/common/back.gif
(版主的素数计算程序没有给出计算孪生素数, 希望提供此功能)
我的素数计算程序是我开发HugeCalc过程中的“副产品”，
由于发现网上有不少人苦于没有足够大的“素数表”而犯愁，
所以才决定把它独立出来与大家分享。
相对来说，孪生素数没有素数的适用面广，用户的需求不是很大（但部分用户可能非常迫切），
所以就没有继续向这方面走，当然还有从程序的整体架构、输入输出界面等方面综合考量。
楼主提出“擂台”，最好是有个参比对象，比如一段代码、或一个程序等，
这样大家才能有目标，并在追逐极限过程中享受快乐，同时达到相互促进。 就我所知，直接产生孪生素数表并没有特别好的办法，可能其复杂度本来就同产生素数表相同。这可能也是为何通常我们看不到专门介绍产生孪生素数表的方法
回复 3# 的帖子
//欢迎大家测试// PI2 = 6810670
// AMD 3600+ 2.0G1.6 s
// intel PD 9403.2G 1.8 s
/*****************************************************************
describition:
input two number beg, end (0 < beg <= end < 2^32)
claculate number of twin primes in the interval
******************************************************************/
# include <stdio.h>
# include <time.h>
# include <math.h>
# include <memory.h>
# include <assert.h>
# include <stdlib.h>
# 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 MASK7
# else
typedef ushort utype;
# define COMP 5
#define MASK15
# 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
utype bits2;
//num of twin primes
int Prime;
//primes in interval
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;
//the table the fist SP primes is removed
/*****************************************************************/
char tablep30;
char tabletp30;
char multp2;
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;
#ifdef _WIN32
# include <windows.h>
DWORD WINAPI Win32ThreadFunc(LPVOID pinfo)
#else
# include <pthread.h>
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.beg = beg;
int tsize = (end - beg) / threadnums;
for (i = 1; i < threadnums; i++)
Threadparam.beg = Threadparam.end = Threadparam.beg + tsize;
Threadparam.end = end;
#ifdef _WIN32
HANDLE thredHand;
DWORD threadID;
for (i = 0; i < threadnums; i++){
thredHand = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE)Win32ThreadFunc, (LPVOID)(&Threadparam), 0, &threadID);
if (thredHand == NULL)
printf("create Win32 thread error = %d\n", GetLastError());
}
WaitForMultipleObjects(threadnums, thredHand, true, INFINITE);
for (i = 0; i < threadnums; i++){
primenums += Threadparam.primenums;
CloseHandle(thredHand);
}
#else
pthread_t pid;
for (i = 0; i < threadnums; i++){
int error = pthread_create(&pid, NULL, POSIXThreadFunc, &Threadparam);
if ( error != 0 )
printf("create pthread error = %d\n", error);
}
for (i = 0; i < threadnums; i++){
pthread_join(pid, NULL);
primenums += Threadparam.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); ){
if (tabletp30 > 0 && pgap != lastgap){
multp2 = pgap - lastgap;
lastgap = pgap;
}
pgap += 2;
mod30 = (twinp30 + pgap * pri30) % 30;
}
}
void inline initMask(utype mask[], int start, uint len)
{
int srid = start % BLOCKSIZE;
memcpy(mask, BaseTpl, sizeof(BaseTpl));
if (start < 32)
*((ushort *)&mask) = ~0xcb6e;
if (srid)
mask |= (MASKN(srid) - 1);
if (len != BLOCKSIZE)
mask |= ~(MASKN(len) - 1);
}
int piRange2(int start, uint len = BLOCKSIZE)
{
static char prime2;
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;
initMask(mask, start--, len);
for (int i = SP + 1, beg = srid - 1, p = Prime; p < sqrtn; p = Prime[++i]){
if (ok){
next = beg + p - start % p;
if ((next & 1) == 0)
next += p;
}
else
next = p * p;
int buf;
const uchar nexti[] = { 1, 2, 3, 4, 5, 0};
while ( tabletp30 == 0 )
next += p * 2;
int st = 5, ps = tablep30 - 1 + (tabletp30 - 1) * 8;
for (int k = 0; k < 6; k++)
buf = p * multp2;
for (; next < len; next += buf])
mask |= MASKN(next);
}
int size = len >> COMP;
for (int k = (srid >> COMP); k <= size; k++){
primenums += bits2];
if (mask < (1 << (TBITSIZE - 1)) && !(mask & 1))
primenums++;
}
int blocki = ++start / BLOCKSIZE;
if (!srid && !(mask & 1)){
prime2 |= 1;
if (prime2 > 1 /*&& blocki > 0*/)
primenums++;
}
if (len == BLOCKSIZE && !(mask[(BLOCKSIZE >> COMP)] & MASKN(BLOCKSIZE - 1)))
prime2 |= 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;
if (pcache == 0)
pcache = piRange2(i * BLOCKSIZE);
primenums += pcache;
}
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++;
mask2 <<= 1;
}
}
//init tablep30
for (i = 0; i < 8; i++)
tablep30 % 30] = i + 1;
//init twinp30
for (i = 0; i < 6; i++)
tabletp30 % 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);
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 : 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 = 2;
for (p = 3; p < maxp; p += 2){
if ( !(BaseTpl & MASKN(p)) ){
Prime[++primes] = p;
for (uint j = p * p; j < maxp; j += p << 1)
BaseTpl |= MASKN(j);
}
}
memset(BaseTpl, 0, sizeof(BaseTpl));
for (int i = 2; i <= SP; i++){
for (int j = Prime, p = j; j < BLOCKSIZE; j += p << 1)
BaseTpl |= MASKN(j);
}
Prime = maxn; Prime = primes;
BaseTpl |= ~(MASKN(BLOCKSIZE) - 1);
BaseTpl[(BLOCKSIZE >> COMP) + 1] = ~0;
return primes;
}
to tprime：
刚才发现你的帖子后半部是斜体，估计是“[i]”被误认为成Discuz!代码了，我给你编辑了下，只要发帖时在“禁用 Discuz!代码”前打钩即可。 原来你是只需要计算数目，而不需要找出所有的素数对 找出来实现也很容易, 位数组已经筛掉不满足条件的，考虑到需要较多内存就没实现.
回复 7# 的帖子
原帖有一些小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)不知道是否有比筛选法更好的？ 我见过的计算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，如果你愿意，在你的程序定型后，可以考虑加入到数学研发网原创软件中进行统一宣传，因为我的收费主页空间还有些剩余（注意：因空间限制，只会收录一些优秀的、精巧的原创作品）。
回复 9# 的帖子
再次更新的算法, 在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)
页:
[1]
2