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

初期値








反復式









πの算出



64bit整数型を用いて円周率を計算する

64bit整数型で小数点を表現する

Visual C++の64bit整数型はUINT64です。
表現できる数値の範囲は、0~18,446,744,073,709,551,615(10進数で20桁弱)です。
C++での掛け算の結果64bitを超えないようにするには、非乗数・乗数ともに32bit(0~4,294,967,295(10進数で10桁弱))に収まっていれば64bitを超えるとこはありません。
上記の式での円周率の計算では、整数部は1桁もあればよいので、小数点を整数に換算して計算するために1,000,000,000を掛けて計算します。
64bitアセンブラの場合 CPU命令で、64bit*64bit=128bit、128bit/64bit=64bitの計算ができるので、64bitの整数の桁数を有効に活用できます。

平方根の計算

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

初期値は

計算例

0 ,  a0= 1000000000 , b0=707106781 , t0=250000000 , p0=1
1 ,  a1= 853553390 , b1=840896415 , t1=228553391 , p1=2
2 ,  a2= 847224902 , b2=847201266 , t2=228473292 , p2=4
3 ,  a3= 847213084 , b3=847213083 , t3=228473292 , p3=8
4 ,  a4= 847213083 , b4=847213083 , t4=228473292 , p4=16
5 ,  a5= 847213083 , b5=847213083 , t5=228473292 , p5=32
3141592619

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

ソースコード

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

#include 
#include 
#include 

#define INT_OFFSET  (UINT64(1000000000))

UINT64 sqrt64(UINT64 src);

int _tmain(void){
	UINT64 a = INT_OFFSET;
	UINT64 a0=0;
	UINT64 b;
	UINT64 t=INT_OFFSET/4;
	UINT64 p=1;
	UINT64 pi;

	b = sqrt64( 2 * INT_OFFSET * INT_OFFSET );
	b = INT_OFFSET*INT_OFFSET / b;

	int n=0;
	_tprintf(_TEXT("%i ,  a%i= %lld , b%i=%lld , t%i=%lld , p%i=%lld\n"),n,n,a,n,b,n,t,n,p);

	do{
		a0 = a;
		a = ( a + b ) / 2;
		b = sqrt64( a0 * b );
		t = t - ( p * (a0 - a) * (a0 - a) ) / INT_OFFSET;
		p = 2 * p;
		++n;
		_tprintf(_TEXT("%i ,  a%i= %lld , b%i=%lld , t%i=%lld , p%i=%lld\n"),n,n,a,n,b,n,t,n,p);
	}while(a0 != a);

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

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