计算孪生素数对

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 <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 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 <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[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)