primenumber..cpp

一个求素数个数的程序 算pi(1000000000)不要1秒每个块如果10 0000, 需要 < 100000次全搜索加上8000 * 8000 < 6400 0000 总的次数小于700

#include <iostream.h>
#include <math.h>
#include <windows.h>

const unsigned default_n=4000000000;
const int maxlen=1000;

int primes[maxlen];
int len;
int sum;
unsigned int n;
int m=7;
int Q=1;
int phiQ=1;
int* v;

bool string2int(unsigned int& n,char* s);
void init(void);
unsigned int phi(unsigned int x,int a);

void main(void){	
	char number_string[80];
	unsigned int time;

	cout<<"calc pi(n)\n";
	cout<<"n(default = "<<default_n<<") = ";
	cin.get(number_string,80);
	if(!string2int(n,number_string))n=default_n;
	
	time=GetTickCount();
	init();
	int num=(int)phi(n,len)-sum+len-1;
	time=GetTickCount()-time;
	GlobalFree(v);
	cout<<"pi(n)="<<num<<"("<<time<<"ms)\n";
	cout<<"press any key to continue...";
	cin.get();
	cin.get();
}

void init(void){
	int max;
	int sqrt_max;
	bool* mark;
	int i,j;
	int len2,len3;
	int s;

	max=(int)(pow(1.0*n,2.0/3.0)+1.5);
	sqrt_max=(int)(sqrt(max)+0.5);
	mark=(bool*)GlobalAlloc(GMEM_FIXED|GMEM_ZEROINIT,max);
	
	for(i=2;i<sqrt_max;++i){
		if(!mark[i]){
			j=i*i;
			while(j<max){
				mark[j]=true;
				j+=i;
			}
		}
	}

	for(len=0,i=2;(unsigned)i*i*i<=n;++i){
		if(!mark[i])primes[++len]=i;
	}

	for(j=max-1,sum=0,s=0,len2=len;(unsigned)i*i<=n;++i){
		if(!mark[i]){
			++len2;
			while((unsigned)i*j>n){
				if(!mark[j])++s;
				--j;
			}
			sum+=s;
		}
	}

	for(len3=len2;i<max;++i){
		if(!mark[i])++len3;
	}

	sum=(len2-len)*len3-sum;
	sum+=len*(len-1)/2-len2*(len2-1)/2;

	if(m>len)m=len;
	for(i=1;i<=m;++i){
		Q*=primes[i];
		phiQ*=primes[i]-1;
	}

	v=(int*)GlobalAlloc(GMEM_FIXED,Q*sizeof(int));
	for(i=0;i<Q;++i)v[i]=i;
	for(i=1;i<=m;++i){
		for(j=Q-1;j>=0;--j){
			v[j]-=v[j/primes[i]];
		}
	}
	
	GlobalFree(mark);
}

bool string2int(unsigned int& n,char* s){
	unsigned int val=0;
	for(int i=0;s[i];++i){
		if(s[i]>'9'||s[i]<'0')return false;
		val*=10;
		val+=s[i]-'0';
		if(val>default_n)return false;
	}
	n=val;
	return true;
}

unsigned int phi(unsigned int x,int a){
	if(a==m){
		return (x/Q)*phiQ+v[x%Q];
	}
	if(x<(unsigned)primes[a]){
		return 1;
	}
	return phi(x,a-1)-phi(x/primes[a],a-1);
}