多倍長で円周率を計算する

■多倍長で円周率を計算する

円周率を計算する公式は沢山ありますが

Machinの公式を使って求めたいと思います

公式で使われているアークタンジェントを求めるには、

■まずは上の公式を使い浮動小数点型を使って円周率を求めるプログラムを書いてみます


#include <stdio.h>

double arctan(double x){
	double a=x;//乗数
	double b=1.0;//分母
	int flag=0;//引き算/足し算を交互に行うフラグ
	double pos;
	double ans=x;

	do{
		a=a*x*x;
		b+=2;
		pos=a/b;

		if(flag=(!flag)){
			ans-=pos;
		}else{
			ans+=pos;
		}
	}while(pos>0.00001);//小数点以下5桁の精度であきらめる
	return ans;
}


int main(){
	printf("%.13lf\n", 4.0*(4.0 * arctan(1.0/5.0) - arctan( 1.0 / 239.0)));
	return 0;
}

▼実行結果

$ ./test
3.1415917721832

浮動小数点型変数を使い円周率が求められました
しかし問題発生、"10進数で大きな桁数の計算を行う"で作成した多倍長の計算方法では小数点型の計算ができません
そこで、分数の分子を大きな数にして解決します

1/3=0.333333333

の計算だと、小数点数が扱えないため計算できません
そのため

100000/3=033333

このようにすると、小数点以下の計算が浮動小数点型を使わなくても出来ることになります
これを固定小数点と言うらしい。

次に計算式を見てみます

どうみても分子が分数でしかも、n乗とかになっています
このままの式では固定小数点での計算方法が思いつきません
いろいろ考えたあげく、式を変形してみました

この式なら分子が正数ですので計算できそうです
そこで、上で作成した浮動小数点型のプログラムを整数型に修正してみます


#include <stdio.h>

int test(int x){
	int a=x;//乗数
	int b=1;//分母
	int flag=0;//引き算/足し算を交互に行うフラグ
	int pos;
	int ans;
	
	ans=10000000/x;
	do{
		a=a*x*x;
		b+=2;
		pos=10000000/(b*a);

		if(flag=(!flag)){
			ans-=pos;
		}else{
			ans+=pos;
		}
	}while(pos>1);

	return ans;
}


int main(){
	printf("%d", 4*(4 * test(5) - test(239)));
	return 0;
}

▼実行結果

$ ./test
31415932

浮動小数点を使わずに円周率が計算できました
これを"10進数で大きな桁数の計算を行う"で作成した多倍長の計算方法で計算するように書き換えるのみです。

■多倍長の計算方法のメイン関数部分を上のプログラムと同じように修正します


#define SIZE 304 //600桁求める

short	a[SIZE],x[SIZE]//乗数
		,b[SIZE]//分母
		,c[SIZE]//分子
		,ans[SIZE]
		,ans2[SIZE]
		,pos[SIZE]
		,pos2[SIZE]
		,val_2[SIZE]={2}// 2
		,val_4[SIZE]={4};// 4

int main(){
	int i;
	int flag=0;//引き算/足し算を交互に行うフラグ

	b[0]=1;//分母
	c[SIZE-1]=1;//分子
	x[0]=5;//乗数
	for(i=0;i<SIZE;i++) a[i]=x[i];
	div(c,x,ans,SIZE);
	for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定
	
	do{
		mul(a,x,pos2,SIZE);
		mul(x,pos2,a,SIZE);
		add(b,val_2,b,SIZE);
		mul(b,a,pos2,SIZE);

		div(c,pos2,pos,SIZE);
		for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

		if(flag=(!flag)){
			sub(ans,pos,ans,SIZE);
		}else{
			add(ans,pos,ans,SIZE);
		}
	}while(cmp0(pos,SIZE)!=0);


	flag=0;//引き算/足し算を交互に行うフラグ

	for(i=0;i<SIZE;i++){ c[i]=0;x[i]=0;b[i]=0;}
	b[0]=1;//分母
	c[SIZE-1]=1;//分子
	x[0]=39;//乗数
	x[1]=2;
	for(i=0;i<SIZE;i++) a[i]=x[i];
	div(c,x,ans2,SIZE);
	for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定
	
	do{
		mul(a,x,pos2,SIZE);
		mul(x,pos2,a,SIZE);
		add(b,val_2,b,SIZE);
		mul(b,a,pos2,SIZE);

		div(c,pos2,pos,SIZE);
		for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

		if(flag=(!flag)){
			sub(ans2,pos,ans2,SIZE);
		}else{
			add(ans2,pos,ans2,SIZE);
		}
	}while(cmp0(pos,SIZE)!=0);

	mul(ans,val_4,pos,SIZE);
	sub(pos,ans2,pos,SIZE);
	mul(pos,val_4,pos2,SIZE);


	dp_print(pos2,SIZE,3);//誤差が生じるため下6桁を切り捨てて固定小数点形式で表示
	printf("\n");

	return 0;
}

▼実行結果

3.
1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196
4428810975 6659334461 2847564823 3786783165 2712019091
4564856692 3460348610 4543266482 1339360726 0249141273
7245870066 0631558817 4881520920 9628292540 9171536436
7892590360 0113305305 4882046652 1384146951 9415116094
3305727036 5759591953 0921861173 8193261179 3105118548
0744623799 6274956735 1885752724 8912279381 8301194912
9833673362 4406566430 8602139494 6395224737 1907021798
6094370277 0539217176 2931767523 8467481846 7669405132

計算速度を考慮していないため、処理速度は激しく遅いですが計算できました

600桁の計算で処理時間が5秒もかかります

PentiumM 1.4GHz メモリ500M　コンパイラVC6　WinXP

■円周率600桁を計算する総てのソースコード


#include <stdio.h>

//  ar+ar2=ar3
void add(const short*ar,const short*ar2,short*ar3,const int size);
//  ar-ar2=ar3
void sub(const short*ar,short*ar2,short*ar3,const int size);
//  ar*ar2=ar3
void mul(const short*ar,const short*ar2,short*ar3,int size);
//  ar/ar2=ar3...ar
void div(short*ar,short*ar2,short*ar3,const int size);
//  printf(ar)
void ar_print(short*ar,const int size);
//配列を固定小数点型として画面表示する ar[size-1].ar[size-2]~ar[low]
void dp_print(short*ar,const int size,const int low);

void turn(short*ar,const int size);
int cmp(const short*ar,const short*ar1,const int size);
int cmp0(const short*ar,const int size);
void r_shift(short*ar,const int size);
void l_shift(short*ar,const int size);




#define SIZE 304 //600桁求める

short	a[SIZE],x[SIZE]//乗数
		,b[SIZE]//分母
		,c[SIZE]//分子
		,ans[SIZE]
		,ans2[SIZE]
		,pos[SIZE]
		,pos2[SIZE]
		,val_2[SIZE]={2}// 2
		,val_4[SIZE]={4};// 4

int main(){
	int i;
	int flag=0;//引き算/足し算を交互に行うフラグ

	b[0]=1;//分母
	c[SIZE-1]=1;//分子
	x[0]=5;//乗数
	for(i=0;i<SIZE;i++) a[i]=x[i];
	div(c,x,ans,SIZE);
	for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定
	
	do{
		mul(a,x,pos2,SIZE);
		mul(x,pos2,a,SIZE);
		add(b,val_2,b,SIZE);
		mul(b,a,pos2,SIZE);

		div(c,pos2,pos,SIZE);
		for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

		if(flag=(!flag)){
			sub(ans,pos,ans,SIZE);
		}else{
			add(ans,pos,ans,SIZE);
		}
	}while(cmp0(pos,SIZE)!=0);


	flag=0;//引き算/足し算を交互に行うフラグ

	for(i=0;i<SIZE;i++){ c[i]=0;x[i]=0;b[i]=0;}
	b[0]=1;//分母
	c[SIZE-1]=1;//分子
	x[0]=39;//乗数
	x[1]=2;
	for(i=0;i<SIZE;i++) a[i]=x[i];
	div(c,x,ans2,SIZE);
	for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定
	
	do{
		mul(a,x,pos2,SIZE);
		mul(x,pos2,a,SIZE);
		add(b,val_2,b,SIZE);
		mul(b,a,pos2,SIZE);

		div(c,pos2,pos,SIZE);
		for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

		if(flag=(!flag)){
			sub(ans2,pos,ans2,SIZE);
		}else{
			add(ans2,pos,ans2,SIZE);
		}
	}while(cmp0(pos,SIZE)!=0);

	mul(ans,val_4,pos,SIZE);
	sub(pos,ans2,pos,SIZE);
	mul(pos,val_4,pos2,SIZE);


	dp_print(pos2,SIZE,3);//誤差が生じるため下6桁を切り捨てて固定小数点形式で表示
	printf("\n");

	return 0;
}

//ar+ar2=ar3
void add(const short*ar,const short*ar2,short*ar3,const int size){
	int i,p_flag;
	for(i=0;i<size;i++){
		ar3[i]=ar[i]+ar2[i];
	}
	for(i=0;i<size;i++){
		if(ar3[i]>99){
			ar3[i]-=100;
			ar3[i+1]+=1;
		}else if(ar3[i]<-99){
			ar3[i]-=-100;
			ar3[i+1]+=-1;
		}
	}

	p_flag=cmp0(ar3,size);
	if(p_flag==0) return;

	for(i=0;i<size;i++){
		if(p_flag>0){
			if(0>ar3[i]){
				ar3[i+1]-=1;
				ar3[i]+=100;
			}
		}else{
			if(0<ar3[i]){
				ar3[i+1]-=-1;
				ar3[i]+=-100;
			}
		}
	}
}

//ar-ar2=ar3
void sub(const short*ar,short*ar2,short*ar3,const int size){
	turn(ar2,size);
	add(ar,ar2,ar3,size);
	turn(ar2,size);
}

//ar*ar2=ar3
void mul(const short*ar,const short*ar2,short*ar3,int size){
	int i,j,pos;
	for(i=0;i<size;i++){
		ar3[i]=0;
	}
	for(i=0;i<size;i++){
		if(ar2[i]){
			for(j=0;j<size;j++){
				ar3[j+i]+=ar[j]*ar2[i];
				//3桁目を繰り上げる
				if(ar3[j+i]>99 || ar3[j+i]<-99){
					pos=ar3[j+i]/100;	
					ar3[j+i]-=(pos*100);
					ar3[j+1+i]+=pos;
				}
			}
		}
	}
}

//  ar/ar2=ar3...ar
//  arが書き換えられて余りになる
void div(short*ar,short*ar2,short*ar3,const int size){
	int i,pos=0;
	int div_flag;//被除数の符号フラグ
	int met_flag;//法の符号フラグ

	for(i=0;i<size;i++){
		ar3[i]=0;
	}

	div_flag=cmp0(ar,size);
	met_flag=cmp0(ar2,size);

	if(!div_flag || !met_flag) return;//法、被除数のどちらかが0であるため計算中止
	if(-1==div_flag) turn(ar,size);//被除数が負数であるため正数にする
	if(-1==met_flag) turn(ar2,size);//法が負数であるため正数にする

	while(1){
		if(cmp(ar,ar2,size)==1 && (ar2[size-1]/10)==0){
			pos++;
			l_shift(ar2,size);
		}else break;
	}

	while(pos>=0){
		if(cmp(ar,ar2,size)>=0){
			sub(ar,ar2,ar,size);
			if(pos%2){
				ar3[pos/2]+=10;
			}else{
				ar3[pos/2]++;
			}
		}else{
			if(pos) r_shift(ar2,size);
			pos--;
		}
	}

	if(-1==div_flag) turn(ar,size);//被除数の符号を戻す
	if(-1==met_flag) turn(ar2,size);//法の符号を戻す
	if(div_flag!=met_flag) turn(ar3,size);//被除数と法の符号が違うため、答えをマイナスにする
}

//画面表示
void ar_print(short*ar,const int size){
	int i;
	int start_flag=0;
	int p_flag=cmp0(ar,size);
	if(p_flag==0){
		printf("0");
		return;
	}
	if(p_flag==-1){
		printf("-");
		turn(ar,size);
	}
	for(i=size-1;i>=0;i--){
		if(ar[i] || start_flag){
			if(!start_flag) printf("%d",ar[i]);
			else printf("%02d",ar[i]);
			start_flag=1;
		}
	}
	if(p_flag==-1) turn(ar,size);
}

//配列を固定小数点型として画面表示する ar[size-1].ar[size-2]~ar[low]
void dp_print(short*ar,const int size,const int low){
	int i,j=0,k=0;
	int p_flag=cmp0(ar,size);
	if(p_flag==0){
		printf("0");
		return;
	}
	if(p_flag==-1){
		printf("-");
		turn(ar,size);
	}

	printf("%d.",ar[size-1]);

	for(i=size-2;i>=low;i--){
		if(j++%5==0) printf(" ");
		if(k++%25==0) printf("\n");
		printf("%02d",ar[i]);

	}
	if(p_flag==-1) turn(ar,size);
}

//符号の反転
void turn(short*ar,const int size){
	int i;
	for(i=0;i<size;i++){
		ar[i]=-ar[i];
	}
}

//ar<ar1 の場合は-1を返す
//ar>ar1 の場合は+1を返す
//ar==ar1 の場合は0を返す
int cmp(const short*ar,const short*ar1,const int size){
	int i;
	for(i=size-1;i>=0;i--){
		if(ar[i]<ar1[i]) return -1;
		if(ar[i]>ar1[i]) return 1;
	}
	return 0;
}

//arがプラスの値をとる場合は+1を返す
//arがマイナスの値をとる場合は-1を返す
//arがゼロの値をとる場合は0を返す
int cmp0(const short*ar,const int size){
	int i;
	for(i=size-1;i>=0;i--){
		if(ar[i]<0) return -1;
		if(ar[i]>0) return 1;
	}
	return 0;
}


//右1シフト
void r_shift(short*ar,const int size){
	int i,pos;
	ar[0]=ar[0]/10;
	for(i=1;i<size;i++){
		pos=ar[i];
		ar[i]=ar[i]/10;
		ar[i-1]+=((pos-(ar[i]*10))*10);
	}

}

//左1シフト
void l_shift(short*ar,const int size){
	int i,pos;
	for(i=size-1;i>=0;i--){
		ar[i]=ar[i]*10;
		pos=ar[i]/100;
		ar[i]-=(pos*100);
		if(pos) ar[i+1]+=pos;
	}

}

▲トップページ > プログラミングの実験