多倍長計算ライブラリGMPで円周率を求める。

多倍長計算したいなー、と思ってググってみるとGMPなるライブラリがあったので、試してみる。

The GNU MP Bignum Library


触り方も何も分からんかったので、この辺を参照。

http://h2np.net/pi/


どうやら

mpf_t 多倍長変数(構造体?)を宣言し
mpf_init で変数の初期化をして


以下の関数で
mpf_add 第一引数に第二引数と第三引数の和を
mpf_sub 同様に差を
mpf_mul 同様に積を
mpf_div 同様に商を代入する


らしい。


どうせなので参考ページとは変えてガウスルジャンドルアルゴリズムで組んでみる。

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

#include 
#include 

#define LOG_TEN_TWO  3.32192809488736234789
#define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2)

int main(int argc,char *argv[])
{
    mpf_t a,b,t,x,y,p,t1,t2,t3,t4;
    long i; 

    long int prec, dprec;
    dprec = 1000000L;
    prec = bprec(dprec);
    mpf_set_default_prec (prec);

    mpf_init(a);
    mpf_init(b);
    mpf_init(t);
    mpf_init(x);
    mpf_init(y);
    mpf_init(p);
    mpf_init(t1);
    mpf_init(t2);
    mpf_init(t3);
    mpf_init(t4);

    mpf_set_ui(a,1); //a = 1

    mpf_set_ui(t1,1); 
    mpf_div_ui(t2,t1,2);
    mpf_sqrt(b,t2); //b = sqrt(1/2)

    mpf_div_ui(t,t1,4); //t = 1/4

    mpf_set_ui(p,1); //p = 1

    for(i = 0;i<10;i++)
    {
        mpf_add(t1,a,b);
        mpf_div_ui(x,t1,2); //x = (a+b)/2

        mpf_mul(t1,a,b);
        mpf_sqrt(y,t1); //y = sqrt(a*b)

        mpf_sub(t1,a,x);
        mpf_pow_ui(t2,t1,2);
        mpf_mul(t3,p,t2);
        mpf_sub(t4,t,t3);
        mpf_set(t,t4); //t = t-p(a-x)^2
        
        mpf_set(a,x); //a = x
        mpf_set(b,y); //b = y

        mpf_mul_ui(t1,p,2);
        mpf_set(p,t1); //p = 2*p

    }

    mpf_add(t1,a,b);
    mpf_mul(t2,t1,t1);
    mpf_mul_ui(t3,t,4);
    mpf_div(t4,t2,t3);

    mpf_out_str(stdout,10,1000000,t4);
    printf("\n");

    return 0;
}


調子に乗って100万桁とかにしているが、ループ回数が少なすぎて、そこまで精度が出ていない気がするorz
(出力はしっかり100万桁出しているが。あとでちゃんとトレーススクリプト書いてみようっと。・・・2786桁しかあってなかったorzもっとループを増やせば!)
とりあえずGMPの概要が分かったので、ここまで。


makefile

CC = gcc
GMP = /usr/lib/libgmp.a

pi:pi.c $(GMP)
    $(CC) -o $@ $< $(GMP) 


近頃Winマシンでもviのキーバインド癖が出るまでになってきました。
マウスジェスチャー初体験以来のイラつきですねw