GMPの紹介

何となくπを多数桁（100万桁）のオーダーで計算するのに、今のPCだとどれくらいかかるんだろうと思い、ネットをググっていたら、アルゴリズムやプログラムの紹介をしているサイトがありました。
そこで使われているライブラリがGMP(Gnu Multipule Precision Arithmetic Library)でして、以前少し耳にはしていましたが詳しいことを調べてなかったので、ちょっと見てみました。
（注：πの多数桁は計算プログラムがサイトに乗っており、その方の計算時間、自分で試した計算時間を見ても1億桁やっても10分ちょっとでできてしまうことがわかり、興味がなくなりました。ただメモリの消費量を見ていると、32bit OSではギリギリに近くなってます。これ以上の桁数の計算をしようと思うと、64bit OSにしないと。）

さてGMPですが、任意の精度で数学計算を行うライブラリです。（もちろんメモリが許す限り、という大前提はつきますが。）　https://gmplib.org/ に本家があります。最新が2013.9更新になっていますから、まだ活動は続いているようです。
とりあえずインストールの説明をします。

１．ソースをDLしてくる。
私は、gmp-5.1.3.tar.bz2を落としてきました。リリースノートを見てみると、CPU（ここではintel, ARMが対象になっていました）毎に最適化をしてより高速化を目標にしているようです。

２．tarボールの解凍、ビルド、インストール
適当な作業用ディレクトリを作り、そこで先ほどDLしてきたtarボールを解凍します。今回、自分はPCとして、Ubuntu 13.10、Mac OSX Mavericksの２台で試してみました。（Mac OSXの方はインストールに少し難があると本家サイトの中に記載がありましたが、特に問題はおきませんでした。ただ困ったのは、OSXでは/usr/local以下にインストールしたものをデフォルトではコンパイル対象にしてくれなくてサンプルプログラムのコンパイルに少し手間がかかりました。）

$ ./configure --enable-cxx　　（c++でも使えるようにオプション追加です）

$ make

$ make check　　　　　　　　（単体試験を多数やってくれます）

$ sudo make install　　　　　（/usr/local以下にライブラリをインストールします）

３．テストプログラム
雰囲気が味わえる程度の簡単なものを示します。

（１）整数の加算
まず基本として整数の加算です。

#include 
#include 

int main (int argc, char* argv[]) {

  mpz_t a, b, c;
  mpz_init(a);
  mpz_init(b);
  mpz_init(c);

  mpz_set_ui(a, 12345);
  mpz_set_str(b, "12345678910987654321", 10);

  mpz_add(c, a, b);
  mpz_out_str(stdout, 10, a);
  printf("\n");

  mpz_out_str(stdout, 10, b);
  printf("\n");

  mpz_out_str(stdout, 10, c);
  printf("\n");

  mpz_clear(a);
  mpz_clear(b);
  mpz_clear(c);

  return 0;
}

最初に整数変数a,b,cを定義しています。（特にprecision（精度）は規定していませんが、整数の場合は必要なだけ確保されるようです）
面白いのが、"mpz_set_str()"関数ですが、通常のC言語では扱えないような多数桁の整数を初期値に使えるように、文字列で設定できる関数です。c=a+bを行い、最後にa,b,cを出力しています。

・実行結果
12345
12345678910987654321
12345678910987666666


ちゃんと加算されていることがわかります。"mpz_out_str()"が出力関数なんですが、第１引数がFILE*で今は標準出力にしていますが、もちろんファイルへの出力にも使えます。また第２引数の"10"ですが、これは基数をいくつで出力するかです。（といっても、8,10,16くらいしか使わないと思いますが）

（２）階乗の計算
多数桁の計算が自由にできるということで、サンプルで階乗の計算プログラムを示します。

#include 
#include 

int main (int argc, char* argv[]) {

  unsigned long int i;
  mpz_t a;

  mpz_init(a);
  mpz_set_ui(a, 1);

  for (i = 1; i <= 100; i++) {
    mpz_mul_ui(a, a, i);
  }

  printf("%ld! = ", (i-1));
  mpz_out_str(stdout, 10, a);
  printf("\n");

  mpz_clear(a);
  return 0;
}

・実行結果
100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

今回は結果の表示を考えて100までにしましたが、10000でも簡単に計算してくれます。

（３）浮動小数点の計算
基本は整数と変わりません。ただGMPの内部では浮動小数点は指数部は通常の浮動小数点形式と同じ（つまり固定長）、その代わり仮数部が任意の桁数にできます。また、変数を宣言／初期化する時に、その仮数部の桁数（正確にはbit長）を指定してやらないといけません。

・例
int digits1 = 20;
int prec1 = digits1 * log2(10);

mpf_t a;

mpf_set_default_prec(prec1);
mpf_init(a);


例に示したものでは変数aの仮数部を２０桁で持つように設定してます。（mpf_set_default_prec()で設定するときbit長で指定しないといけないので、２０桁が何ビットになるか計算してます）

但し、いくら仮数部を多数桁もてたとしても、指数部の極端に異なる浮動小数点同士で演算を行えば桁落ちは発生します。そのあたりの注意事項は変わりません。（仮数部をかなり大きく持てるので、計算誤差は大分減らせるでしょうが）

（４）分数の計算
GMPではなんと分数の計算もできます。

#include 
#include 

int main (int argc, char* argv[]) {
 mpq_t a, b, c;

 mpq_init(a);
 mpq_init(b);
 mpq_init(c);

 mpq_set_si(a, 4, 7);
 mpq_set_si(b, 3, 8);

 mpq_add(c, a, b);

 mpq_out_str(stdout, 10, a);
 printf("\n");

 mpq_out_str(stdout, 10, b);
 printf("\n");

 mpq_out_str(stdout, 10, c);
 printf("\n");

 return 0;
}

内部的には整数の分子と分母で持っているのかな。上記は、4/7 + 3/8の計算をしています。

・実行結果
4/7
3/8
53/56


こんな感じで、ちゃんと分数形式で出力もしてくれます。当然、分数で計算後、浮動小数点への変換をしてくれる関数も持っていますので、やたら沢山の（整数の）分数を掛けたり、割ったりしなければいけない計算式の計算がかなり誤差を気にしなくてもよくなります。

さてこんな便利なGMPライブラリなんですが、残念ながら完全にはWindowsに対応していません。独自にVisual Studioのprojを作っている人がいましたが、make check（単体試験）まで対応できているんだろうか？Windowsにインストールしている人はほとんどcygwinでした。（Linux, Mac使っている人間が今更cygwinを使う必要性が殆ど無いので、これは試していません。）
configureをCMakeあたりにしてくれると面白いんですが、それでも単体試験の問題は残るでしょうし。（本家のサイト読んでたら、とにかくmake checkは必ずやれとうるさく書いてありました。OSの環境、CPUのバージョンの微妙な違いでFPUの結果が変わったりしますから、そこでライブラリの保証をしないと、危なくて使えないモノになってしまいますしね。）