prime_bit.cpp

素数统计快速算法

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <memory.h>

typedef unsigned int uint;

#ifndef N5
	#define BLOCKSIZE 150150 
	// 2*3*5*7*11*13 * 5 = 150150
	#define ST 6
#else
	#define BLOCKSIZE 510510 
// 2*3*5*7*11*13 * 17=510510
	#define ST 7
#endif

#define MOVE 5

// 需要优化
// 1。 汇编编写部分执行较多的代码
// 2。 改为多线程
// 3。 更新算法， 速度达到提高10 - 100 倍。

int main()
{
	int i, j, p;
	//	unsigned int QRT; //the square root of the end of current block
	uint Start; //the first multiple of a prime in current block
	clock_t tstart;
	uint BasePrime[BLOCKSIZE / 10] = { 2 }; //store the first prime block
	unsigned __int64 PrimeCnt = 1; //prime counter
	unsigned __int64 KPoint = 0;   //the start of current block
	uint BaseTpl[(BLOCKSIZE >> MOVE) + 2]; //the block template
	uint CurBuf[(BLOCKSIZE  >> MOVE) + 2]; //current block
	uint umask[32];
	for (i = 0; i < 32; i++)
		umask[i] = 1u << i;
	memset(BaseTpl, -1u, sizeof(BaseTpl));

	BaseTpl[0] = ~umask[1];
	for (i = 1; i * i < BLOCKSIZE; i += 2){
		if (BaseTpl[i >> MOVE] & umask[i & 31] ){
			BasePrime[PrimeCnt++] = i;
			for (j = i * i; j < BLOCKSIZE; j += i * 2)
				BaseTpl[j >> MOVE] &= ~umask[j & 31];
		}
	}

	while (i < BLOCKSIZE){
		if (BaseTpl[i >> MOVE] & umask[i & 31])
			BasePrime[PrimeCnt++] = i;
		i += 2;
	}
	//	printf("prime number = %I64d, max primer = %d \n",PrimeCnt,BasePrime[PrimeCnt-1] );
	memset(BaseTpl, -1u, sizeof(BaseTpl));
	for (i = 1; i < ST; i++){
		for (p = BasePrime[i]; p < BLOCKSIZE; p += BasePrime[i]* 2)
			BaseTpl[p >> MOVE] &= ~umask[p & 31];
	}

	uint LQRT = (uint)sqrt(BLOCKSIZE);

	for (j = 1; j < BLOCKSIZE; j++){
		memcpy(CurBuf, BaseTpl, sizeof(CurBuf));
		KPoint += BLOCKSIZE;
		const uint QRT = (uint)sqrt(KPoint + BLOCKSIZE) + 1;
		if (j % 1000 == 0)
			printf(">%04dms\t计算到%11I64u\t查出素数:%10I64u个\n", clock() - tstart , KPoint, PrimeCnt);

		for (i = ST, p = BasePrime[i]; p < QRT; p = BasePrime[++i]){ 
			if (p < LQRT)
				Start = (uint)(p - (KPoint - 1) % p) - 1;
			else
				Start = (unsigned __int64)p * p - KPoint;
			if ((Start & 0x1) == 0)
				Start += p;
			while (Start < BLOCKSIZE){
				CurBuf[Start >> MOVE] &= ~umask[Start & 31];
				Start += p << 1;
			}
		}
		LQRT = KPoint;
		for (i = 1; i < BLOCKSIZE; i += 2){
			if (CurBuf[i >> MOVE] & umask[i & 31])
				//if(PrimeCnt < BLOCKSIZE) BasePrime[PrimeCnt] = KPoint + i;
				PrimeCnt++;
		}
	}
	return 0;
}