p2(big).cpp

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

////////////////////////////////////////////////////////////////////////////////////////////////
//	P2.cpp																					  //
//  作用：用于计算函数P2																	  //
//	编写者：初旭（00512038）																  //
//  完成时间：2006、5、24																	  //
////////////////////////////////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "conio.h"

__int64 p2big(int* PList, __int64 x,int x_13,int a)
{
	__int64* Block = (__int64*)malloc(sizeof(__int64)*x_13);				// 用于存储每个区间
	__int64* Prime_in_Block = (__int64*)malloc(sizeof(__int64)*x_13);		// 用于存储区间中的素数
	__int64* Pi_of_Block = (__int64*)malloc(sizeof(__int64)*(x_13+1));		// 用于存储区间中每个数的PI值
	__int64* noteP = (__int64*)malloc(sizeof(__int64)*x_13);				// 用于筛选记录

	__int64 Delta;
	__int64 S;
	__int64 Previous_Pi=0;
	__int64 result=0;
	
	__int64 x_12 = pow(x,0.5);							// x的1/2次方
	__int64 x_16 = pow(x_13,0.5);						// x的1/6次方
	__int64 x_14 = pow(x,0.25);							// x的1/4次方
	__int64 note = x_12 / x_13;							// 记录在哪个区间出现x的1/2次方
	__int64	PIx_12;										// x的1/2次方的PI值
	__int64 i,j,k,t,n,l,r;
	__int64 length;										// 不整区间的长度
	double u= pow(x,2.0/3);
	__int64 x_23 = (__int64)(u + 0.000001);					// x的2/3次方
	__int64 Bl;											// 区间的第一个数
	__int64 block;										// 区间的最后一个数
	__int64 I0;

	if((x_12 % x_13) != 0)
		note++;
	
	for(t=1 ; t<=x_13+1 ; t++)
	{
/////////////////////////////////////////////////////////////////////////////////////////////////
		if(t==1)					// 另行处理（这时已经有素数表了）
		{
			for (i =0; i<x_13; i++)
				Block[i] = i+1;
			noteP[0] = x_13/2+1;
			for (i = 1; i<a; i++)
				noteP[i] = (x_13-PList[i])/(PList[i]<<1) +1;	// 这里是为筛选做准备工作

			Previous_Pi = a;
			S=0;
		}

//////////////////////////////////////////////////////////////////////////////////////////////////
//	下面的这个区间情形与正常区间的处理相同，注释略												//
//////////////////////////////////////////////////////////////////////////////////////////////////
		else if(t==x_13 +1)							// 这个区间另行处理
		{
			length = x_23 - (t-1)*x_13;
			if (length <= 0)
				continue;

			free (Block);
			free (Prime_in_Block);
			free (Pi_of_Block);
				
			__int64* Block = (__int64*)malloc(sizeof(__int64)*length);
			__int64* Prime_in_Block = (__int64*)malloc(sizeof(__int64)*length);
			__int64* Pi_of_Block = (__int64*)malloc(sizeof(__int64)*(length+1));

			for (i =0; i<length; i++)
				Block[i] =(t-1)*x_13 + i +1;
			Bl= (t-1)*x_13;
			block = Block[length -1];

			// 1.用p<=x_13 的素数筛选这个区间,筛到的置为0 ///////////////////////////
			
			for(k= noteP[0]<<1; k <= block; k+=2)			// 用2正常的筛选一遍
			{
				noteP[0]++;									// 每筛一次，倍数加1
				Block[k -Bl -1] = 0;
			}
	
			for (i=1; i<a; i++)								// 其他素数用2倍筛选
			{
				for(k= ((noteP[i]*PList[i])<<1)+PList[i]; k <= block; k+=(PList[i]<<1))
				{
					noteP[i]++;
					Block[k -Bl -1] = 0;
				}
			}


			// 2.给Prime_in_Block 赋值////////////////////////////////////////////////
			j=0;
			for (i=0; i<length; i++)
			{
				if(Block[i] !=0)
				{
					Prime_in_Block[j] = Block[i];
					j++;
				}
				else
					Block[i]=i + (t-1)*x_13 +1;
			}
			__int64 NumPrime_in_Block = j;
				
			// 3.通过Prime_in_Block 和pi((i-1)N)来计算：Pi_of_Block,
			//	 也就是{Pi(y):(i-1)N <=y <=iN}									///
			j = 0;
			n = 0;
			for (i=0; i<NumPrime_in_Block;)
			{
				for(;( Block[j] <= Prime_in_Block[i] || i==NumPrime_in_Block) && j<length; j++)
				{
					if(Block[j] == Prime_in_Block[i])
					{
						n++;
						Pi_of_Block[j]= n+Previous_Pi;
						j++;
						i++;
					}
					Pi_of_Block[j] = n+Previous_Pi;
				}
			}

			// 4. 定位区间I[i] 中的所有素数（I[i] 就是原文中的I[j],I[j]的范围与j的关系见原文）//
			if(x_13 > ( k = (x/(t*x_13 + 1))))
				l=x_13;
			else
				l=k;

			if((k= (x/((t-1)*x_13+1))) < x_12)
				r=k;
			else
				r=x_12;

			if( l>=r )
				continue;

			__int64* I= (__int64*)malloc(sizeof(__int64)*(r-l+1));

			for(i=0; i <(r-l); i++)
			{
				I[i] = l + i + 1;
			}
			I0=I[0];

			for(i=0; PList[i]<=x_14; i++)
			{
				for(k=(l/PList[i]+1)*PList[i]; k<=r; k+=PList[i])
					I[k-I0]=0;
			}

			// 5.计算表{ y = [x/p] : p在I[i]中（上面定位的）},通过与Pi_in_Block比较，
			//   求出x/p_in_Block表中的y 的pi 值（可以直接写入这个表中将原来的值覆盖），
			//   将x/p_in_Block中的元素求和，就是Delta
			Delta = 0;
			for(i=0; i<(r-l); i++)
			{
				if (I[i] == 0 )
					continue;
				else
				{
					I[i] = x / I[i];
					I[i] = Pi_of_Block[I[i] -(t-1)*x_13 -1 ];
					Delta +=(__int64)(I[i]);
				}
			}
			free(I);
			S+=Delta;

			free (Block);
			free (Prime_in_Block);
			free (Pi_of_Block);
		}
/////////////////////////////////////////////////////////////////////////////////////////////////
/////////////				以上为最后一个区间的 处理，以下为一般区间的处理//////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////

		else						// 一般的情况
		{
			for (i =0; i<x_13; i++)
				Block[i] +=x_13;							// 给区间赋值
			
			Bl= (t-1)*x_13;
			block = Block[x_13 -1];

			// 1.用p<=x_13 的素数筛选这个区间,筛到的置为0 ///////////////////////////
			for(k= noteP[0]<<1; k <= block; k+=2)			// 用2正常的筛选一遍
			{
				noteP[0]++;									// 每筛一次，倍数加1
				Block[k -Bl -1] = 0;
			}
	
			for (i=1; i<a; i++)								// 其他素数用2倍筛选
			{
				for(k= ((noteP[i]*PList[i])<<1)+PList[i]; k <= block; k+=(PList[i]<<1))
				{
					noteP[i]++;
					Block[k -Bl -1] = 0;
				}
			}

			// 2.给Prime_in_Block 赋值////////////////////////////////////////////////
			j=0;
			for (i=0; i<x_13; i++)
			{
				if(Block[i] !=0)
				{
					Prime_in_Block[j] = Block[i];			// 给这个区间的素数数组赋值
					j++;
				}
				else
					Block[i]=i + (t-1)*x_13 +1;
			}
			__int64 NumPrime_in_Block = j;

			// 3.通过Prime_in_Block 和pi((i-1)N)来计算：Pi_of_Block,			///
			//	 也就是{Pi(y):(i-1)N <=y <=iN}									///
			j = 0;
			n = 0;

			for (i=0; i<NumPrime_in_Block;)
			{
				for(;( Block[j] <= Prime_in_Block[i] || i==NumPrime_in_Block) && j<x_13; j++)
				{	
					if(Block[j] == Prime_in_Block[i])		// 如果区间里的这个数值是素数
					{
						n++;								// 那么从此开始PI值加1
						Pi_of_Block[j]= n+Previous_Pi;		
						j++;
						i++;
					}
	
					Pi_of_Block[j] = n+Previous_Pi;
				}
			}
			Previous_Pi = Pi_of_Block[x_13 - 1];

			if(t == note)									// 如果这个区间包含x的1/2次方
			{
				PIx_12= Pi_of_Block[x_12 - (t-1)*x_13 -1];	// 那么记录x的1/2次方的PI值
			}

			// 4. 定位区间I[i] 中的所有素数（I[i] 就是原文中的I[j],I[j]的范围与j的关系见原文）//
			if(x_13 > (k =x/(t*x_13 + 1)))
				l=x_13;
			else
				l=k;

			if((k= (x/((t-1)*x_13+1))) < x_12)
				r=k;
			else
				r=x_12;

			if( l>=r )										// 以上是求出I区间的上下限
				continue;

			__int64* I= (__int64*)malloc(sizeof(__int64)*(r-l));

			for(i=0; i <(r-l); i++)
			{
				I[i] = l + i + 1;							// I区间赋值
			}

			I0=I[0];
			for(i=0; PList[i]<=x_14 && i<a; i++)			// 在I区间中筛选素数
			{
				for(k=(l/PList[i]+1)*PList[i]; k<=r; k+=PList[i])
					I[k-I0]=0;
			}


			// 5.计算表{ y = [x/p] : p在I[i]中（上面定位的）},通过与Pi_in_Block比较，
			//   求出x/p_in_Block表中的y 的pi 值（可以直接写入这个表中将原来的值覆盖），
			//   将x/p_in_Block中的元素求和，就是Delta
			Delta = 0;
			for(i=0; i<(r-l); i++)
			{
				if (I[i] == 0 )
					continue;
				else
				{
					I[i] = x / I[i];
					I[i] = Pi_of_Block[I[i] - (t-1)*x_13 -1 ];	// 对照求出PI值
					Delta +=(__int64)(I[i]);
				}
			}
			free(I);

			S+=Delta;
		}
	}

	result = a*(a-1)/2 - ((__int64)PIx_12)*(PIx_12 - 1)/2 +S;

	return (result);

}