super pi を作ろう

円周率πを100万桁まで計算するプログラムを作ってみました。super pi もどきです。ソースコードと解説を載せます。

ソースコードは80行程度の非常に短いものですが、このような簡単な実装でも最近のPCなら100万桁の円周率が一分も掛らず計算できます。

ポイント

要点は、

  • 計算精度が問題となるので、GMP(任意精度計算ライブラリ)を使う
  • 計算アルゴリズムは、Square AGM *1と呼ばれるものを使う

です。

GMPについては、http://d.hatena.ne.jp/pyopyopyo/20090303 をご覧下さい。

ソースコード

/**
 * pi.cpp
 *
 * written by pyopyopyo at gmail dot com
 */
#include <math.h>
#include <gmpxx.h>
#include <iostream>

// N桁の円周率を計算
#define  N  100*100*100

using namespace std;

int
main()
{
     // N桁の数値計算に必要な精度を prec(bits) とすると
     //
     // prec = log2(10^N) 
     //      = N * log2(10)

     // 計算誤差を考えるのがめんどうなので,とりあえず(N+5)桁で計算してみる
     int prec = (N+5)*log2(10);
     mpf_set_default_prec(prec);

     fprintf(stderr, "prec:= %d bits\n", prec);

     const mpf_class two ("2.0");

     mpf_class a ("1.0");
     mpf_class b (a/sqrt(mpf_class("2.0")));
     mpf_class t (a);
     mpf_class x (a);


     mpf_div_ui(t.get_mpf_t(), t.get_mpf_t(), 4);

     mpf_class tmp(a);

     const int P = ceil(log2(N));
     for (int k=1; k<=P; ++k) {
          tmp = a;
          a = a + b;
          mpf_div_ui(a.get_mpf_t(), a.get_mpf_t(), 2);

          b = sqrt( tmp*b );

          tmp -= a;
          mpf_pow_ui(tmp.get_mpf_t(), tmp.get_mpf_t(), 2);
          t -= x * tmp;

          x *= two;
     }

     tmp = a + b;
     mpf_pow_ui(tmp.get_mpf_t(), tmp.get_mpf_t(), 2);
     tmp /= t;
     mpf_div_ui(tmp.get_mpf_t(), tmp.get_mpf_t(), 4);


     mp_exp_t expo;
     char *str = mpf_get_str(0, &expo, 10, N+1, tmp.get_mpf_t());

     for (int i=expo,count=1; i<N; i+=10, ++count) {
          char tmp[12]={0};
          strncpy(tmp, &str[i], 10);
          cout << tmp ;
          if (0 == (count%10))
               cout << endl;
          else
               cout <<" ";
     }

     return 0;
}

コンパイル方法と、実行方法

$ gcc -lgmpxx pi.cpp -o pi
$ time ./pi

計算時間の目安としては

ぐらいになります。

なお出力は、以下のURLで公開されている100万桁の円周率πのデータの形式に合わせています。

ftp://pi.super-computing.org/pub/pi10m/pi10m.ascii.01of10

解説

  • 計算精度の指定
    int prec = (N+5)*log2(10);
    mpf_set_default_prec(prec);

gmpのデフォルトの計算精度を、(N+5)桁に設定しています。

log2(10^(N+5)) → (N+5)*log2(10) ということです。

なお、今回はN=100万なので、 precは約332万になります。つまり、GMPの変数は、一変数につき 400Kbyte 程度のメモリを消費することになります。

  • #include

gmpc++インターフェイスを利用しています。

たとえば、gmpのcインタフェースで書いた以下のコードは

 mpf_t a, b, c;
 mpf_init(a);
 mpf_init(b);
 mpf_init(c);

 mpf_set_str(a, "1234", 10);
 mpf_set_str(b, "5678", 10);
 mpf_add(c, a, b);

 mpf_out_str(stdout, 10, 10, c);

c++インタフェースでは

 mpf_class a,b,c;
 a="1234";
 b="5678";
 c=a+b;
 cout << c;

と簡潔に記述できます。

さらに

 mpf_class a("1234"),b("5678"),c;
 c=a+b;
 cout << c;

とも書け、さらに圧縮して

 cout << mpf_class("1234")+mpf_class("5678");

と書くことも出来ます。


数値計算の常識
数値計算の常識
posted with amazlet at 09.03.09
伊理 正夫 藤野 和建
共立出版
売り上げランキング: 112092
おすすめ度の平均: 5.0
5 意外と忘れがちだったり、重要だけど知らなかったりするのが 常識なんです。
5 文字どおり「常識」にしたい内容
5 困る前に

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ
William H. Press William T. Vetterling Saul A. Teukolsky Brian P. Flannery
技術評論社
売り上げランキング: 9168
おすすめ度の平均: 4.0
4 教育的価値>実用的価値
5 最高です
4 第二版の翻訳がほしい
4 他の数値計算の参考書を買う前に
4 アルゴリズムの宝庫!

πの歴史 (ちくま学芸文庫)
ペートル ベックマン
筑摩書房
売り上げランキング: 152704