s2.cpp

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

#include "stdafx.h"




pF initF(bint x_13,int *PList)
{
	//动态开辟结构体F的数组	arrayF[x_13]
	pFarrayF;
	arrayF=(pF)malloc(sizeof(F)*x_13);
	bint i,j;//??j可能为int
	int* p;
	for(i=0;i<x_13;i++)
	{
		arrayF[i].m=i;
		
		arrayF[i].miu_m=//??在其他函数中已被计算出来
		
		for(p=PList,j=0;j<a;j++)
		{
			if(i%*(p+j-1)==0)
			{
				arrayF[i].f_m=*(p+j-1);
				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;
}


bint s2(bint x,bint a,int *PList)
{
	bint s2=0;//函数返回值
	
	//初始化F
	pF arrayF;
	arrayF=F(bint x_13);
	//-----------------------------------------------3.7

	bint N=x_13;//x_13表示x的1/3次方取整
	//据此求k的变化范围

	//动态建立fai[a]  //a=pi(x_13)

	//-----------------------------------------------3.8

	//动态建立a[floor(log2N)][floor(N/2^i)+1]  //??可以采取链表数组形式创建a

	//a[floor(log2N)][floor(N/2^i)+1]=2^i      //初始化数组a
	//-----------------------------------------------3.9

	
	int* p=PList; //传入p表
	
	bint k,/*搜索区间参数*/j/*fai的第二个参数*/;

	bint y,l;//l从1到N变化//??正常区间，非正常区间待定

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

	bint m;//父结点
	
	
	for(k=1;k<=N+2;k++)     //搜索第k个区间
	{
		for(j=1;j<=a;j++)	//fai的第二个参数的变化
		{
		//1.固定k,j的一次筛选，更新数组a
			while(*p<=*(p+j-1))
			{
				for(l=1;l<=N;l++)//??当Bk是正常区间，非正常区间待定
				{
					y=(k-1)*N+l;
					if(y%(*p)==0)
					{
						//所有I_i,j包含y的a[i][j]--
					}
				}

				p++;
			}
			//--------------------------------------------------------3.9.5


		//2.估算叶子值
			L=max((x/((k*N+1)*(*(p+j)))),1);
			U=min((x/((((k-1)*N+1)*(*(p+j)))),N);
		    //---------------------------------------------------------3.13

			for(m=L;m<=U;m++)
			{
				if(arrayF[m].f_m>*(p+j)&&arrayF[m].miu_m!=0)
				{
					y=floor(x/(m*(*(p+j))));
					l=y-(k-1)*N;
					//提取l的二进制表示，即提出一个数组e[m-1],??可以考虑用位操作
					bint fai_yj=0,h;
					for(i=0;i<m;i++)
					{
						for(j=0,h=1;j<i;j++)
						{
							h+=pow(2,e[j]-e[i]);
						}
						fai_yj+=a[e[i]][h];
					}
					
					fai_yj+=fai[j];
			//---------------------------------------------------------3.10
					s2+=miu_m*fai_yj;//s2的数值累积
				}
			}
			


		//3.更新fai表
			l=N;
			//提取l的二进制表示，即提出一个数组e[m-1]
			bint fai_yj=0,h;
			for(i=0;i<m;i++)
			{
				for(j=0,h=1;j<i;j++)
				{
					h+=pow(2,e[j]-e[i]);
				}
				fai_yj+=a[e[i]][h];
			}
					
			fai_yj+=fai[j];
			fai[j]=fai_yj;//3.8的更新
			//---------------------------------------------------------3.10

		}
	}
	

	
	
	
	return s2;
}