gb-bbp.cpp

bbp算法计算圆周率精确到千万位以上

//BBP法计算圆周率
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define H	16

//定义全局变量
int binary[33];				//整形数组存储二进制数
long double dpi[1000];		//存储圆周率，每个单元存储8位16进制小数
long double	time=0;

//计算余数
double res(double);

//计算第n项的8位16进制小数
double hexcal(int);

//将8位16进制小数取为整的函数，并防止其防止越界
double wint(double);


main()
{
	int num,hex,dec,n,i,j,temp1,temp2;
	double shex=0,dhex,htemp;
    long int dtemp;
	
	//输入要计算的最大小数位
	printf("\n***********************BBP法计算圆周率*****************************\n");
    do
	{
		printf("\n输入要计算的最大小数位：\n");
		scanf("%d",&num);
		if(num<=0)
			printf("输入错误！请输入一个非负整数！\n");
	} while(num<=0);

	dhex=pow(16,8);			//dhex赋值为16的8次方
	
	for(i=0;i<1000;i++)		//存储圆周率的数组初始化为0
		dpi[i]=0;
	
	hex=int(num/9.63)+1;		//十进制位n换成十六进制位n2
	
	for(i=1;i<=hex;i++)		
	{
		n=8*(i-1);
		shex=hexcal(n);			//依次计算第i项的8位16进制小数
		htemp=wint(shex*dhex);	//将其转化为整数，取8位（十六进制类型的）
		dpi[i]=htemp/dhex;
	}

	dec=(int(1.20411998*8*hex)+4)/5;		//转化为十进制输出
	printf("\n\n*************************计算结果显示*****************************\n\n");
	for(j=1;j<=dec;j++)                     
	{   	
		dpi[hex]=100000*dpi[hex];
		for(i=hex;i>1;i--)
		{
			dtemp=long(dpi[i]);
			dpi[i]=dpi[i]-dtemp;
			dpi[i-1]=100000*dpi[i-1]+dtemp/dhex;
			time++;
		}
		dtemp=long(dpi[1]);
		dpi[1]=dpi[1]-dtemp;
		if((j-1)%10==0)
		{	
			temp1=(j-1)*5+1;
			temp2=temp1+49;
			if((j-1)==0)
				printf("\t pi=3.");
			else 
				printf("\n%5d-%5d位:",temp1,temp2);
		}
		if(dtemp<10)
			printf(" 0000%d",dtemp);
		else 
			if(dtemp<100)
				printf(" 000%d",dtemp);
			else
				if(dtemp<1000)
					printf(" 00%d",dtemp);
				else		
					if(dtemp<10000)
						printf(" 0%d",dtemp);
					else
						printf(" %d",dtemp);
	}

	printf("\n\n*****************************计算复杂度********************************\n\n");
	printf("时间复杂度:%12.f\n",time);
	printf("空间复杂度:      1000(long double)\n");

	return 0;
}

double res(double p)
{
	double wh=16,wtemp=1;
	int j;
	for(j=1;j<binary[0]+1;j++)
	{
		if(binary[j]==1)
		{	
			wtemp=wh*wtemp;
			wtemp=wtemp-p*int(wtemp/p);		//取模的余数
		}
		if(j==binary[0]) break;
		wh=wh*wh;
		wh=wh-p*int(wh/p);					//取模的余数
	}
    
	return wtemp/p;							//返回余数
}

double hexcal(int n)
{
	int i,j,ntemp;
	long double hex0=1,hex2,hex3,hex4,hex16=1,thex,sh=0;

	for(i=0;i<n+13;i++)
	{
		hex2=hex0+3;
		hex3=hex0+4;
		hex4=hex0+5;
		if(i<n+1)
		{   
			if(i==0)
				sh=2.0/15.0;
			else 
			{
				ntemp=n-i;
				for(j=1;j<33;j++)   //把m-k变为二进制数
				{
					time++;
					binary[j]=ntemp%2;
					ntemp=ntemp/2;
					if(ntemp==0) break;
				}
				binary[0]=j;
				sh=sh+(4*res(hex0)-2*res(hex2)-res(hex3)-res(hex4));//计算所有小数和
				while(sh>1)  sh=sh-1;			//保证s为纯小数
				while(sh<0)	 sh=sh+1;
			}
		}
		else
		{
			time++;
			hex16=hex16/16;					//计算没有整数部分的后12项
			thex=hex16*(4.0/hex0-2.0/hex2-1.0/hex3-1.0/hex4);
			sh=sh+thex;
			while(sh>1) sh=sh-1;
			while(sh<0) sh=sh+1;			
			if(thex<1e-17)  break;
		}
		hex0=hex0+8;
	}
	return sh;
}

double wint(double n)
{
	double sh,sh1,sh2;
	long th1,th2;
	sh=n/100000.0;
	th1=int(sh);
	sh1=sh-th1;
	sh2=sh1*100000;
	th2=int(sh2);

	return	th1*100000.0+(long)sh2;
}