ガウス=ルジャンドルのアルゴリズム

初期値








反復式









πの算出



double型を用いて円周率を計算する

浮動小数点 double型

Visual C++のdouble型は-1.79769313486231570E+308~1.79769313486231570E+308を表せます。
仮数部は10進数で約18桁表現できます。INT64とdouble型はメモリ上では8バイトで同一であるが、浮動小数点の強みで精度の維持が可能である

平方根の計算

Aの平方根は下記の式で反復すると近似できる

初期値は

計算例

0 , a0=1.0000000000000000 , b0=0.7071067811865476 , t0=0.2500000000000000 , p0=  1
1 , a1=0.8535533905932737 , b1=0.8408964152537146 , t1=0.2285533905932738 , p1=  2
2 , a2=0.8472249029234942 , b2=0.8472012667468916 , t2=0.2284732910809006 , p2=  4
3 , a3=0.8472130848351929 , b3=0.8472130847527655 , t3=0.2284732905222318 , p3=  8
4 , a4=0.8472130847939792 , b4=0.8472130847939792 , t4=0.2284732905222318 , p4= 16
5 , a5=0.8472130847939792 , b5=0.8472130847939792 , t5=0.2284732905222318 , p5= 32
3.1415926535897940

小数点以下14桁目まで正確な値が計算されている。

ソースコード

// ガウス=ルジャンドルのアルゴリズムによるdouble型を使用した円周率の算定

#include 
#include 
#include 
#include 

#define INT_OFFSET  (double(1))

double sqrtf(double src);

int _tmain(void){
	double a = 1;
	double a0 = 0;
	double b;
	double t=1.0/4.0;
	double p=1;
	double pi;
	double temp;
	b = sqrtf( 2.0 );
	b = 1 / b;

	int n=0;

	do{
		a0 = a;
		a = ( a + b ) / 2.0;
		b = sqrtf( a0 * b );
		t = t - ( temp = p * (a0 - a) * (a0 - a) );
		p = 2.0 * p;
		_tprintf(_TEXT("%i , a%i=%18.16f , b%i=%18.16f , t%i=%18.16f , p%i=%3.0f\n"),n,n,a,n,b,n,t,n,p);
		++n;
	}while(a0 != a);

	pi=((a+b))*((a+b))/(4*t);
	_tprintf(_TEXT("%18.16f\n"),pi);
	return 0;
}

double sqrtf(double src){
	double ans=src/2;
	double ans0;
	double t;
	do{
		ans0=ans;
		t=src/ans0;
		ans=(t+ans0)/2;
	}while(ans0 != ans);
	return ans;
}