s2.cpp

用Meisell-Lehmer算法计算PI(x)

#include "stdafx.h"

typedef __int64 bint;

#include "s2.h"

#include "stdlib.h"	

#include "stdio.h"

#include "time.h"

#include "math.h"

int a;


	
pF initF(int x_13,int *PList,int *Miu)    //对结构体数组初始化的函数
{
	
	pF arrayF;              //动态开辟结构体F的数组	arrayF[x_13]

	arrayF=(pF)malloc(sizeof(F)*x_13);   //指针类型为pF

	int i,j;
	int* p;
	for(i=0;i<x_13;i++)
	{
		arrayF[i].m=i+1;           //给m赋值
		
		arrayF[i].miu_m=Miu[i];    //m相应的miu值
		
		for(p=PList,j=0;j<a;j++)   //该循环用于找m的最小质因数
		{
			if(i==0)
			{
				arrayF[i].f_m=0;
			}

			else if((i+1)%p[j]==0)
			{
				arrayF[i].f_m=p[j];

				break;
			}
		}
	}

	
	return arrayF;

}


bint max(bint x,bint y)
{
	if(x>=y)
		return x;
	else
		return y;
}
bint min(bint x,bint y)
{
	if(x>=y)
		return y;
	else
		return x;
}


int * e_sz(int l,int * m)      //提取l的二进制形式并将其二进制指数存在e_sz数组中，m同时记录数组元素个数
{
	int *e1;

	e1=(int*)malloc(sizeof(int)*32);     //整型的l最多有32位，故用长度为32的数组

	unsigned int a=2147483648;          //位运算时用的模版数

	int i,j=0;

	for(i=sizeof(int)*8-1,*m=0;i>=0;i--)
	{
		if(l&a)
		{
			e1[j]=i;
			j++;
			(*m)++;
		}
		a=a>>1;
	}
	return e1;					

}			

bint s2(bint x,int a,int *PList,int x_13,int *Miu)
{
	bint s2=0;    //函数返回值

	pF arrayF;     //算法中的F表

	arrayF=initF(x_13,PList,Miu);

	int N3;     //不整区间的长度

	double n2 = pow(x,2.0/3);
	bint N2 = (bint)(n2+0.0000001);        //精确求得x的三分之二次方取整
	bint N=x_13;       //x_13表示x的1/3次方取整

	int i,j1;//循环控制变量

	bint * fai;        
	      
	fai=(bint*)malloc(sizeof(bint)*((bint)(a+1)));
	for(i=0;i<=a;i++)
	{
		fai[i]=0;
	}                     //建立循环过程中需要使用的fai表并初始化
	
	int log2N=floor(log10(N)/log10(2))+1;        //log2N是以后反复要使用的数

	int *Ni=(int*)malloc(sizeof(int)*log2N); 
    
	for(i=0;i<log2N;i++)
	{
		Ni[i]=(N>>i)+2;
	}                               //Ni数组中也是以后反复要使用的数，这样来提速

    int * * a_sz;                   //动态建立算法中的a数组

    a_sz=(int **)malloc(sizeof(int*)*log2N);

	for(i=0;i<log2N;i++)
	{
		a_sz[i]=(int*)malloc(sizeof(int)*Ni[i]);
	}                                           
////////////////////////////////////////////////////////
	
	int*e;
	
	int l, m_e=0;

	bint**lvalue;

	lvalue=(bint**)malloc(sizeof(bint)*N);

	for(l=0;l<N;l++)                //对于每一个从1到N的l值，求得算法中的h的值，并放在lvalue数组中，这样可提速
	{
		e=e_sz(l+1,&m_e);

		lvalue[l]=(bint*)malloc(sizeof(bint)*m_e);

		for(i=0;i<m_e;i++)           //不能写for(i=0,lvalue[l][i]=1;i<N;i++)因lvalue[l][i]=1是变量
		{
			lvalue[l][i]=1;
			for(j1=0;j1<i;j1++)
				lvalue[l][i]+=(1<<(e[j1]-e[i]));
		}
		free(e);
	}

//////////////////////////////////////////////
	int**lj;
	lj=(int**)malloc(sizeof(int)*log2N);
	for(i=0;i<log2N;i++)
	{
		lj[i]=(int*)malloc(sizeof(int)*N);
		for(l=0;l<N;l++)
		{
			lj[i][l]=(l>>i)+1;
		}

	}                          //也是要反复使用的一组数，用于计算a_sz数组第二维的下标

	int *ai=(int*)malloc(sizeof(int)*log2N);   //要反复使用的一组数，用于给a_sz数组赋值

	for(i=0;i<log2N;i++)
	{
		ai[i]=1<<i;
	}

///////////////////////////////////////////////////////////

	int * eN;          //l等于N时的e数组，用于更新fai表
	int m_eN=0;
	eN=e_sz(N,&m_eN);
////////////////////////////////////////////////////////////

	int**el=(int**)malloc(sizeof(int*)*N);       //要反复使用的一组数的一组数，计算父结点时使用，可提速
	int*ml=(int*)malloc(sizeof(int)*N);
	for(l=0;l<N;l++)
	{
		el[l]=e_sz(l+1,&ml[l]);
	}

//////////////////////////////////////////////////////////
	int* p=PList;              //传入p表
	
	bint  k;                   //搜索区间参数

	bint y;                    //实际被筛到的数

	bint L,U;                  //父结点搜索区间的下限与上限

	bint m;                    //父结点
	
	bint j;                    //循环变量
	
	int notel;                 //跳跃筛选参量
	
	int delta;                 //跳跃量

	int * flag;                //表征元素y是否已被筛出过的数组

	flag=(int*)malloc(sizeof(int)*N);

 	int N2_N = N2/N;           //整区间的数目

	for(k=1;k<=N2_N;k++)       //处理第k个区间
	{	
		
//首先寻找j=0时的特殊叶子并计算其值

		L=max(ceil((x/((k*N+0.999999)*2))),1);           
		U=min((x/((((k-1)*N+1.0)*2))),N);					
		for(m=L;m<=U;m++)
		{
			if(arrayF[m-1].miu_m!=0&&arrayF[m-1].f_m>2&&m*2>N)
			{
				y=floor(x/(m*2));
				if(arrayF[m-1].miu_m==-1)
					s2+=y;
				else
					s2-=y;               //将该值累加到s2上
			}
		}              

/////////////////////////////////////////////////////////////////////////////////////////////////

		for(i=0;i<N;i++)                //初始化flag数组
			flag[i]=1;
		for(i=0;i<log2N;i++)            //初始化a_sz数组
		{
			for(j1=0;j1<Ni[i];j1++)
			{
				a_sz[i][j1]=ai[i];
			}
		}
		
//////////////////////////////////////////////////////////////////////////////////////////////////

//j不等于0时的特殊叶子的计算

		for(j=1;j<=a;j++)               //fai的第二个参数的变化
		{	
			l=PList[j-1]-((k-1)*N)%PList[j-1];      //直接计算第一个被筛出的l值

			while(l<=N&&flag[l-1]==0)
			{
				l+=PList[j-1];
			}                            //找到尚未被筛出的l值
			if(j==1)
				delta=2;
			else 
				delta=2*PList[j-1];      //跳跃量

			for(notel=l;notel<=N;notel+=delta)
			{
				if(flag[notel-1])         //当没被筛出时
				{
					for(i=0;i<log2N;i++)
					{
						a_sz[i][((notel-1)>>i)+1]--;
					}
				}
				flag[notel-1]=0;
			}                            //跳跃筛选
	
			L=max(ceil((x/((k*N+0.999999)*PList[j]))),1);
			U=min((x/((((k-1)*N+1.0)*PList[j]))),N);              //求父结点的上下限
		
			for(m=L;m<=U;m++)                           //该循环找父结点计算特殊叶子的值并累加
			{	
				if(j==a)
					break;
				
				if(arrayF[m-1].miu_m!=0 && arrayF[m-1].f_m>PList[j] && m*PList[j]>N)
				{
					y=floor(x/(m*PList[j]));
					l=y-(k-1)*N;
					bint fai_yj=0;   
					for(i=0;i<ml[l-1];i++)
					{
						fai_yj+=a_sz[el[l-1][i]][lvalue[l-1][i]];	
					}
						
					fai_yj+=fai[j];
					if(arrayF[m-1].miu_m==-1)
						s2+=fai_yj;
					else
						s2-=fai_yj;
				}
			}                                   

			bint fai_yj=0;
			for(i=0;i<m_eN;i++)
			{
				fai_yj+=a_sz[eN[i]][lvalue[N-1][i]];
			}
			fai_yj+=fai[j];
			fai[j]=fai_yj;                      //更新fai表
		}
	}
	free(eN);
	free(flag);

	for(i=0;i<log2N;i++)
		free(a_sz[i]);
	free(a_sz);
	
//////////////////////////////////////////////////////////////////////////////////////////


//以下是特殊区间的处理

	if(N3=N2-(k-1)*N)
	{	
		int log2N3=floor((log10(N3)/log10(2)))+1;      //log2N3是以后反复要使用的数

		int *N3i=(int*)malloc(sizeof(int)*log2N3);      //N3i数组中也是以后反复要使用的数，这样来提速

		for(i=0;i<log2N3;i++)
		{
			N3i[i]=(N3>>i)+2;
		}
///////////////////////////////////////////////////////////////////////////////////////////////
		a_sz=(int **)malloc(sizeof(int*)*(log2N3));     //动态建立算法中与该区间相应的a数组

		for(i=0;i<log2N3;i++)
		{	
			a_sz[i]=(int*)malloc(sizeof(int)*N3i[i]);
		}

		flag=(int*)malloc(sizeof(int)*N3);               //建立算法中与该区间相应的flag数组

//////////////////////////////////////////////////////////////////////////////////////////////////

//以下还是先求j=0时的特殊叶子的值并加至s2上

		L=max(ceil((x/((N2+0.999999)*2))),1);
		U=min((x/((((k-1)*N+1.0)*2))),N);					
		for(m=L;m<=U;m++)
		{
			if(arrayF[m-1].miu_m!=0&&arrayF[m-1].f_m>2&&(m*2)>N)
			{
				y=floor(x/(m*2));
				if(arrayF[m-1].miu_m==-1)
					s2+=y;
				else
					s2-=y;
						
			}
		}
//////////////////////////////////////////////////////////////////////////////////////////////////
		for(i=0;i<N3;i++)                  //初始化flag数组
			flag[i]=1;

		for(i=0;i<log2N3;i++)              //初始化a_sz数组
		{
			for(j1=0;j1<N3i[i];j1++)
			{
				a_sz[i][j1]=(1<<i);
			}
		}
//////////////////////////////////////////////////////////////////////////////////////////////////

//j不等于0时的特殊叶子的计算

		for(j=1;j<=a;j++)
		{	
			l=PList[j-1]-((k-1)*N)%PList[j-1];               //直接计算第一个被筛出的l值

			while(l<=N3&&flag[l-1]==0)                       //找到尚未被筛出的l值
			{
				l+=PList[j-1];
			}
			if(l>N3)
				continue;
			if(j==1)                                          //以下与整区间的处理手法同
				delta=2;
			else
				delta=2*PList[j-1];
			for(notel=l;notel<=N3;notel+=delta)
			{
				if(flag[notel-1])
				{
					for(i=0;i<log2N3;i++)
					{
						a_sz[i][((notel-1)>>i)+1]--;
					}
				}
								
				flag[notel-1]=0;
			}

			L=max(ceil((x/((N2+0.999999)*PList[j]))),1);
			U=min((x/((((k-1)*N+1.0)*PList[j]))),N);

			for(m=L;m<=U;m++)
			{
				if(j==a)
					break;
				if(arrayF[m-1].miu_m!=0&&arrayF[m-1].f_m>PList[j]&&m*PList[j]>N)//1的最小质因子怎样？
				{
					y=floor(x/(m*PList[j]));
					l=y-(k-1)*N;
					bint fai_yj=0;
					for(i=0;i<ml[l-1];i++)
					{
						fai_yj+=a_sz[el[l-1][i]][lvalue[l-1][i]];	
					}
					fai_yj+=fai[j];
					if(arrayF[m-1].miu_m==-1)
						s2+=fai_yj;
					else
						s2-=fai_yj;
				}
			}
		}
		free(flag);

		for(i=0;i<log2N3;i++)
			free(a_sz[i]);
		free(a_sz);
	}
///////////////////////////////////////////////

//所有有关数组的free
	
	free(fai);

	free(Ni);

	free(ai);

	for(i=0;i<N;i++)
	{
		free(lvalue[i]);
	}
	free(lvalue);

	for(i=0;i<log2N;i++)
	{
		free(lj[i]);
	}
	free(lj);
/////////////////////////
	for(l=0;l<N;l++)
	{
		free(el[l]);
	}
	free(el);
	free(ml);
//////////////////////////	
	return s2;
}