円周率計算プログラム(64bit整数による計算)(32/64bit)

icon 項目のみ表示/展開表示の切り替え

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

初期値

a_{0}=1

\displaystyle b_{0}=\frac{1}{\sqrt{2}}

\displaystyle t_{0}=\frac{1}{4}

p_{0}=1

反復式

\displaystyle a_{n+1}=\frac{ {a_{n}+b_{n}} }{{2}}

\displaystyle b_{n+1}=\sqrt{ a_{n} b_{n} }

t_{n+1}=t_{n}-p_{n}(a_{n}-a_{n+1})^{2}

p_{n+1}=2p_{n}

πの算出

\pi \approx \frac{ (a+b)^{2} }{4t}

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の平方根は下記の式で反復すると近似できる
a_{n+1}=(A/a_{n}+a_{n})/2
初期値は
a_{0}=A/2

計算例


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 <windows.h>
#include <stdio.h>
#include <tchar.h>

#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;
}