円周率計算プログラム(浮動小数点 double)(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}}

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}

πの算出

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

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

浮動小数点 double型

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

平方根の計算

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

計算例


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

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