孪生素数对php编写,计算孪生素数对 - 编程擂台 - 数学研发论坛 - Powered by Discuz!...

//欢迎大家测试

// 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;

}