プログラムモグモグ

円周率を1億桁計算しました！ ― その試行錯誤の詳しい経緯と結果 ー

雑記

春休み暇ですし, 円周率を計算してみることにしました. エントリーが長くなりましたがお付き合いください.

はじめにお断り

私は円周率計算に関しては全くの素人です. もっとスケーラブルなコードの書き方があると思いますので, あまりここばかりあてにしないでください.

時間がないせっかちな人へ

コード書く試行錯誤をだらだら書いたので, 割とエントリーが長くなっちゃってます. 結論をここに書きます.

• 当初の目標は円周率1000万桁を計算することでした.
• 結局のところは, 後に上げる参考文献を実装しただけです.
• 計算アルゴリズムは, Chudnovskyアルゴリズムです. 最近の円周率計算の記録はこのアルゴリズムに基づいています. Gauss AGMというアルゴリズムを用いている他のソフトウェアと実行時間を比較することで, Chudnovskyアルゴリズムがいかに速いかを示しました.
• 円周率の小数点以下1億桁の計算を行いました. かかった時間は807.20秒, ファイル出力も含めると995.11秒です. 他のサイトの結果と付きあわせて, 100000000桁まで一致することを確認しました.
• [追記:date=2012/03/05]再度計算し, 1億100桁まで計算し, 1億98桁まで正確に求められたことを確認しました. この時の計算時間は816.12秒, ファイル出力も含めると1005.39秒でした.[/追記]

とてもお世話になった文献

Computation of 2700 billion decimal digits of Pi using a Desktop Computer, Fabrice Bellard, Feb 11 2010 (4th revision), http://bellard.org/pi/pi2700e9/pipcrecord.pdf

もし, ご自分で円周率を計算しようと思っておられる方は, 今直ぐこのpdfを開いて, 頑張って実装しましょう. 素晴らしいです. 以下, この記事を｢Bellardさんの記事｣と呼ぶことにします. このエントリーのコードは全てこの記事を見て書いたものです.

And I have to say thank you to Mr.Bellard, for this great article. All the codes and ideas in this entry is based on your article.

パソコンのスペック

三年前に買ったパソコンですので, 皆さんからすればかなり遅いと思います. しかも, マルチスレッド(コア?)とかよく分かんなくて, 結局シングルコアで計算しています.

CPU: Intel(R) Core(TM)2 Duo CPU P8600 @2.40GHz

Memory: 1.8GiB

OS: Ubuntu 11.10

Kernel: Linux 2.0.0-16-generic

皆さんのパソコン上での計算時間のご報告, お待ちしております. このエントリーの最後に掲載したいと思います.

目標を立てる

目標は大事です. モチベーションにつながるからです. 自分は, 最初に次の目標を立てました.

• 円周率1000万桁を計算する.
• それが全て(もしかしたら最後の桁を除いて)正しい円周率と一致することを, どうにかして検算する.
• 高速化する. 100万桁で10秒を切る.
• ブログを書く.

一つ目の1000万桁ですが, テキトーです. 1億はまぁムリダナ(･☓･)...的な感じでした, 私みたいなへっぽこが書けるプログラムじゃぁメモリーが足りなくなりそうだと思っていました. ベタにプログラム書いても, メモリーが余裕で足りそうな感じの桁が1000万桁でした. あと, 日本人には有名な｢スーパーπ｣のサイト*1を見てみますと, あれって3355万桁までしか計算できないんですね. じゃぁ4000万桁行ってみますかな...! みたいな. 結局1億桁行っちゃったんですよね...

二つ目は, どんな手を使ってもいいから, 正しい円周率と比較する. これは本当は, すごく難しい.

三つ目. この目標は, とあるのページのベンチマーク*2を見て立てました. y-cruncher*3などの速いプログラム100万桁で3秒切るんですよね. 標準的なPCの上で.

四つ目. とっても重要. みんなのためよりも, 自分のため.

調べる

とりあえず, 考えていても何もはじまらないので, 調べました. 最初に見たサイトは, *4でした. 5兆桁も計算しているのだからスゴイアルゴリズムに違いない! y-cruncherというプログラムをここで初めて知りました.

y-cruncherはChudnovskyアルゴリズムを用いています*5. で, このアルゴリズムをググります. Wikipedia*6によると,

The Chudnovsky algorithm is a fast method for calculating the digits of π.

http://en.wikipedia.org/wiki/Chudnovsky_algorithm

"a" fast method と, 少し謙遜していますが, 速いんですね! これを使いましょう. また, 2009/12のBellardさん(この方はさっきの記事を書かれた方ですね)と, 2010/8/3のYee，近藤さんの記録はChudnovskyの公式に基づいているらしいです*7. 採用です! 実装しよう!

近藤茂さんは2010年に5兆桁, 2011年に10兆桁を達成された方です*8 *9. いろいろ調べていると, この方のお話がUstreamに上がっていました*10. 私のこのブログを読むよりも, この動画のほうが役に立つ気がします. というかこのポストを見ている人は全員この動画を見に行きなさい!

Chudnovskyアルゴリズム

このアルゴリズムで計算するという方針だけは立ったのですが, 式がなんかややこしいです.

コピペ用: \frac{1}{\pi} = 12 \sum_{n=0}^\infty \frac{(-1)^n (6n)! (A + B n)}{(3n)!(n!)^3 C^{3n + 3/2}} \qquad \left(A = 13591409, B = 545140134 , C = 640320\right)

先ほど挙げたUstreamの動画の近藤茂さんのお話で, このアルゴリズムが, nを1個づつ進めると14桁ずつ精度が上がるということです. よく分からないですが, n=8まで計算すると100桁超えるんですね. スゴイです.

でも, 出てくるのが円周率の逆数ですね. うーん. しかも階乗... まずは手計算してみましょうか. 最初はn=0で.

WolframAlpha用: 640320^(3/2) / 12 / 13591409

これをWolfram Alpha*11で計算してみますと, 3.14159265358973420とかなりました. これは実際の円周率と小数点以下13桁まで一致します!!!

n=1まで求めてみましょう.

WolframAlpha用: 1 / 12 / ( 13591409 / 640320^(3/2) - ( 6! * (13591409 + 545140134) ) / ( 3! * 640320^( 9/2 ) ) )

この値は, 本物の円周率と小数点以下27桁まで一致します*12. 本当に14桁改善しました.

スゴイなぁと思いながら, この公式をぼーっと眺めながら考えていたのですが, ここでBellardさんの記事を見つけました. 素晴らしい記事です. すこしこの記事に基づいて, Chudnovskyアルゴリズムをイジイジしてみましょう.

まずは, の部分を変形します. 階乗はいやだからね.

結果, 階乗割る階乗が, 自然数の掛け算になりました. ここでを定義しましたが, 他に便利な変数を定義しておきます.

大事なのは, これらが整数という事です. AもBも整数なのでは整数, そして, も整数です.

これらを用いて, 最初のChudnovskyの式を変形します.

ここで,

と定義しました.

さて, 式に踊らされずに目標を思い出しましょう. 目標は円周率を計算することです!

ところが, さっきのの定義で行くと, はとてもとても小数です. 効率よく計算するためには, 整数に持ち込まなくてはなりません. じゃぁ, に何かを掛けて整数にすればいいですね. ここで新たにとを次のように定義します.

とをかけ合わせ, それを新たにと定義します.

ももも整数ですから, も整数になりました. このを使うと,

となります.

あとは, , , に関して漸化式を作ります.

初項は

とすれば計算がうまくいきます. もうのことは忘れましょう! 我々はこれらの初項と漸化式から何らかの方法でとを求めれば良いのです. 人によったら最初のが気になるかもしれません. しかし, 実際にコードを書いてみると分かるのですが, ここは最後の最後に掛けるので, 計算時間はなんてことないです.

GMPをインストールする

GMPのサイト*13からソースコードをダウンロードして, インストールしてください.

(実は最初はMPFRを使っていたのですが, この記事を書いてる途中で, もうGMPでいいじゃんって思ったので, 記事を書きながら実装しなおしています...)

次のコードが実行できればOK. ファイル名はtest_gmp.c.

#include <stdio.h>
#include <gmp.h>

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

  mpz_t a, b;

  mpz_init(a);
  mpz_init(b);

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

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

  mpz_clear(a);
  mpz_clear(b);

  return 0;
}

 $ gcc test_gmp.c -lgmp && ./a.out
12345
12345678910987654321

もうちょっと計算してみます.

#include <stdio.h>
#include <gmp.h>

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);
  mpz_out_str(stdout, 10, a);
  printf("\n");

  mpz_clear(a);
  return 0;
}

 $ gcc test_gmp.c -lgmp && ./a.out
100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

うん, 合ってる合ってる*14. 実際書いてみると, C言語のインターフェースでも割と書けるもんですね...

ChudnovskyアルゴリズムをGMPで実装してみる

さっきの漸化式に基づいて書いてみましょう. 後で整数演算に持ち込みますが, まぁ最初ですので全ての変数を浮動小数点(mpf_t)でやってみます. なにも考えずにすみますし, そもそもうまく動くか分かりませんし... GMPに慣れていない方も, ここのコードは眺めれば分かると思います. というか私自身, GMPを使うのは今回が初めてですし, マニュアル*15も分かりやすいですので, 何も困ることはありませんでした.

/* Calculate pi based on Chudnovsky algorithm, using GMP */
/* [1] Computation of 2700 billion decimal digits of Pi using a Desktop Computer,
       Fabrice Bellard, Feb 11 2010 (4th revision),
       http://bellard.org/pi/pi2700e9/pipcrecord.pdf */
#include <stdio.h>
#include <gmp.h>
#include <math.h>

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

  int digits = 100;
  int prec = digits * log2(10);
  int i;
  int n = digits / 14;

  mpf_t pi, p, q, a, P, Q, T, A, B, C;

  /* Initialize Floating point numbers */
  mpf_set_default_prec(prec);
  mpf_init(pi);
  mpf_init(p);
  mpf_init(q);
  mpf_init(a);
  mpf_init(P);
  mpf_init(Q);
  mpf_init(T);
  mpf_init(A);
  mpf_init(B);
  mpf_init(C);

  /* Assignment */
  mpf_set_str(A, "13591409", 10);
  mpf_set_str(B, "545140134", 10);
  mpf_set_str(C, "640320", 10);
  mpf_set_ui(P, 1);
  mpf_set_ui(Q, 1);
  mpf_set_ui(T, 0);

  /* Main loop. about 14 * n digits precision */
  for (i = 1; i <= n; i++) {
    /* p = (2 * i - 1) * (6 * i - 5) * (6 * i - 1) */
    mpf_set_ui(p, (2 * i - 1) * (6 * i - 5) * (6 * i - 1));
    /* q = C^3 * i^3 / 24 */
    mpf_set_ui(q, i * i * i);
    mpf_mul(q, q, C);
    mpf_mul(q, q, C);
    mpf_mul(q, q, C);
    mpf_div_ui(q, q, 24);
    /* a = (-1)^i * (A + B * i) */
    mpf_mul_ui(a, B, i);
    mpf_add(a, a, A);
    if (i & 1) {
      mpf_neg(a, a);
    }
    /* P(0, N) = P(0, N - 1) p */
    mpf_mul(P, P, p);
    /* Q(0, N) = Q(0, N - 1) q */
    mpf_mul(Q, Q, q);
    /* T(0, N) = T(0, N - 1) * q + a * P */
    mpf_mul(T, T, q);
    mpf_mul(a, a, P);
    mpf_add(T, T, a);
  }

  /* pi = C ^ (1 / 2)     */
  mpf_sqrt(pi, C);
  /*      * C             */
  mpf_mul(pi, pi, C);
  /*      * Q             */
  mpf_mul(pi, pi, Q);
  /*      / (T + A * Q)   */
  mpf_mul(Q, Q, A);
  mpf_add(Q, Q, T);
  mpf_div(pi, pi, Q);
  /*      / 12            */
  mpf_div_ui(pi, pi, 12);

  mpf_out_str(stdout, 10, digits, pi);

  mpf_clear(pi);
  mpf_clear(p);
  mpf_clear(q);
  mpf_clear(a);
  mpf_clear(P);
  mpf_clear(Q);
  mpf_clear(T);
  mpf_clear(A);
  mpf_clear(B);
  mpf_clear(C);
  return 0;
}

pi.cとして保存します. 実行すると

 $ gcc pi.c -lgmp && ./a.out
0.3141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068e1

最後に e1 って書いてますので, 多分合ってますね! まずは最初の関門をくぐりました!

評価方法を考える

ちゃんと動いたし, 1000桁, 1万桁ってどんどん行こうぜ! っていう方, ちょっと待って下さいね.

計算結果がどれだけの精度で合っているかを確かめなければなりません. 一番ドン臭いのはソースコードに実際の円周率を書くことですが, さすがにアレですね.

本物の円周率計算の現場では, まだ誰も計算したことのない桁を計算していらっしゃいますので, 実際の円周率の値など分かっていません. ですから別のアルゴリズムで検算するわけですが, ちょっと私には難しそうです.

今回は, 単純に円周率を書いたテキストファイルを用意して, pi.cの出力と引くという方法を取りました. こうしておくと, 他の人が書いたプログラムを実行して評価したり, そういうケースにも便利です. 円周率300桁をpians.txtに書きます.

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273

円周率うん万桁は, 載っているサイトがいくらでもありますので, 適当に拾ってきましょう*16.

引き算をするプログラムを書きます.

/* Substruct two floating point */
/*  $ gcc sub.c -lgmp -o sub */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>

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

  int digits = 1000;
  int prec = digits * log2(10);

  mpf_t x, y;
  FILE* fpx;
  FILE* fpy;

  if (argc != 3) {
    printf("Usage: sub file1 file2");
    exit(EXIT_FAILURE);
  }

  /* Initialize Floating point numbers */
  mpf_set_default_prec(prec);
  mpf_init(x);
  mpf_init(y);

  /* Read from files */
  fpx = fopen(argv[1], "r");
  if (fpx == NULL) {
    printf("Couldn't load %s", argv[1]);
    exit(EXIT_FAILURE);
  }
  fpy = fopen(argv[2], "r");
  if (fpy == NULL) {
    printf("Couldn't load %s", argv[2]);
    exit(EXIT_FAILURE);
  }
  mpf_inp_str(x, fpx, 10);
  mpf_inp_str(y, fpy, 10);

  /* Substruct x = x - y, and output x */
  mpf_sub(x, x, y);
  mpf_out_str(stdout, 10, 20, x);
  printf("\n");

  mpf_clear(x);
  mpf_clear(y);
  return 0;
}

コンパイルして, さっきのpi.cの結果を引きましょう!!!

 $ cat pians.txt
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273
 $ gcc sub.c -lgmp -o sub
 $ gcc pi.c -lgmp && ./a.out > output.txt
 $ ./sub pians.txt output.txt
-0.17851913486717693353e-100

この結果, 円周率を小数点以下100桁の精度で求めることができた, という事が分かりました.

そりゃ当たり前です. さっきのpi.cで

  int digits = 100;

ってしていましたから. ここを500に書きなおし, もう一回計算してみます.

 $ gcc pi.c -lgmp && ./a.out > output.txt
 $ ./sub pians.txt output.txt
-0.72458700660631558817e-300

あれれ??? 500桁までいかないぞ???

この300ってのは, 今度はpians.txtの方ですね. pians.txtを円周率1000桁にして, 再度計算してみます.

 $ cat pians.txt
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
 $ wc pians.txt
   1    1 1001 pians.txt
 $ gcc pi.c -lgmp && ./a.out > output.txt
 $ ./sub pians.txt output.txt
0.29833673362440656643e-499

きちんと500桁まで求まりました. (499桁ですが, こういうのをだいたい500桁とこのエントリーでは言うことにします. digitsを500とすることは, 3.14の3から数えて500桁ということなので, 小数点以下で言うと499桁になってしまいます.)

このようにして考えていくと, 検算に関して次の3つの値に気をつけるべきという事が分かります.

• pians.txt の桁数
• sub.c の digits
• pi.c の digits (と n だが, n はおよそ digits / 14 とすればいい).

計算したい桁までpians.txtとsub.cを書きなおしていけばいいわけです. まぁこれでもどんくさいですがw

pi.cを見直す

さっきのpi.cはあまりにもアレです. すべて浮動小数点演算ですし. それでも何桁まで計算できるかトライしてみると, 5559桁で何故か止まります.

 $ # pians.txt を20000桁にする, sub.c, pi.c のdigitsを20000にし, subをコンパイル
 $ wc pians.txt
   1    1 20001 pians.txt
 $ gcc sub.c -lgmp -o sub
 $ gcc pi.c -lgmp && ./a.out > output.txt
 $ ./sub pians.txt output.txt
-0.35347439649713750729e-5559
 $

怪しいのは次の計算.

    /* p = (2 * i - 1) * (6 * i - 5) * (6 * i - 1) */
    mpf_set_ui(p, (2 * i - 1) * (6 * i - 5) * (6 * i - 1));

溢れてそうです. 書き直します.

    /* p = (2 * i - 1) * (6 * i - 5) * (6 * i - 1) */
    mpf_set_ui(p, (2 * i - 1));
    mpf_mul_ui(p, p, (6 * i - 5));
    mpf_mul_ui(p, p, (6 * i - 1));

 $ ./sub pians.txt output.txt
-0.19999999999999999985e-19999

はい, 20000桁行きました! 2万ですよ2万! だんだん桁数が伸びて行って楽しくなって来ませんか?

しかし, このまま行っても今度は時間がかかりすぎてダメです. 時間を計測するために, ちょっとコードを書き加えます. output.txtへの出力と標準出力を分けます.

#include <time.h>
int main (int argc, char* argv[]) {

  clock_t start, end;
  FILE* fp;

  start = clock();
  // ...
  /* mpf_out_str(stdout, 10, digits, pi); */
  fp = fopen("output.txt", "w");
  if (fp == NULL) {
    printf("couldn't open output.txt");
    exit(EXIT_FAILURE);
  }
  mpf_out_str(fp, 10, digits, pi);

  mpf_clear(pi);
  // ...

  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

  return 0;
}

 $ wc pians.txt
     1      1 200003 pians.txt
 $ gcc sub.c -lgmp -o sub
 $ gcc pi.c -lgmp && ./a.out
21.230 s
 $ ./sub pians.txt output.txt
0.13009369188925586578e-49999
 $

2万桁で2.8秒, 5万桁で20秒です. 三年前のPCでテストしているのですが, 目標は100万桁で10秒, とてもじゃないけどこのままでは行きそうにありません...

整数演算に落とす

実は, , , , , , , , , これらは全て整数です. ですから, ループの中の計算は全て整数で行えるはずです. さっそく書き直しました. mpf_t(浮動小数点)とmpz_t(整数)を繋ぐ関数は, mpf_set_z関数です.

/* Calculate pi based on Chudnovsky algorithm, using GMP */
/* [1] Computation of 2700 billion decimal digits of Pi using a Desktop Computer,
       Fabrice Bellard, Feb 11 2010 (4th revision),
       http://bellard.org/pi/pi2700e9/pipcrecord.pdf */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>

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

  clock_t start, end;
  FILE* fp;

  int digits = 20000;
  int prec = digits * log2(10);
  int i;
  int n = digits / 14;

  start = clock();

  mpf_t pi, temp;
  mpz_t p, q, a, P, Q, T, A, B, C, C3over24;

  /* Initialize GMP numbers */
  mpf_set_default_prec(prec);
  mpf_init(pi);
  mpf_init(temp);
  mpz_init(p);
  mpz_init(q);
  mpz_init(a);
  mpz_init(P);
  mpz_init(Q);
  mpz_init(T);
  mpz_init(A);
  mpz_init(B);
  mpz_init(C);
  mpz_init(C3over24);

  /* Assignment */
  mpz_set_str(A, "13591409", 10);
  mpz_set_str(B, "545140134", 10);
  mpz_set_str(C, "640320", 10);
  mpz_mul(C3over24, C, C);
  mpz_mul(C3over24, C3over24, C);
  mpz_div_ui(C3over24, C3over24, 24);
  mpz_set_ui(P, 1);
  mpz_set_ui(Q, 1);
  mpz_set_ui(T, 0);

  /* Main loop. about 14 * n digits precision */
  for (i = 1; i <= n; i++) {
    /* p = (2 * i - 1) * (6 * i - 5) * (6 * i - 1) */
    mpz_set_ui(p, (2 * i - 1));
    mpz_mul_ui(p, p, (6 * i - 5));
    mpz_mul_ui(p, p, (6 * i - 1));
    /* q = C^3 * i^3 / 24 */
    mpz_set_ui(q, i * i);
    mpz_mul_ui(q, q, i);
    mpz_mul(q, q, C3over24);
    /* a = (-1)^i * (A + B * i) */
    mpz_mul_ui(a, B, i);
    mpz_add(a, a, A);
    if (i & 1) {
      mpz_neg(a, a);
    }
    /* P(0, N) = P(0, N - 1) p */
    mpz_mul(P, P, p);
    /* Q(0, N) = Q(0, N - 1) q */
    mpz_mul(Q, Q, q);
    /* T(0, N) = T(0, N - 1) * q + a * P */
    mpz_mul(T, T, q);
    mpz_mul(a, a, P);
    mpz_add(T, T, a);
  }

  /* pi = C ^ (1 / 2)     */
  mpf_set_z(temp, C);
  mpf_sqrt(pi, temp);
  /*      * C             */
  mpf_mul(pi, pi, temp);
  /*      * Q             */
  mpf_set_z(temp, Q);
  mpf_mul(pi, pi, temp);
  /*      / (T + A * Q)   */
  mpz_mul(Q, A, Q);
  mpz_add(Q, Q, T);
  mpf_set_z(temp, Q);
  mpf_div(pi, pi, temp);
  /*      / 12            */
  mpf_div_ui(pi, pi, 12);

  /* mpf_out_str(stdout, 10, digits, pi); */
  fp = fopen("output.txt", "w");
  if (fp == NULL) {
    printf("couldn't open output.txt");
    exit(EXIT_FAILURE);
  }
  mpf_out_str(fp, 10, digits, pi);

  mpf_clear(pi);
  mpf_clear(temp);
  mpz_clear(p);
  mpz_clear(q);
  mpz_clear(a);
  mpz_clear(P);
  mpz_clear(Q);
  mpz_clear(T);
  mpz_clear(A);
  mpz_clear(B);
  mpz_clear(C);
  mpz_clear(C3over24);

  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

  return 0;
}

実行します.

 $ gcc pi.c -lgmp && ./a.out
0.030 s
 $ ./sub pians.txt output.txt
-0.17961434609008952241e-19999
 $ # digits = 50000に書きなおす
 $ gcc pi.c -lgmp && ./a.out
0.220 s
 $ ./sub pians.txt output.txt
0.13009369188925586578e-49999
 $ # digits = 100000に書きなおす
 $ gcc pi.c -lgmp && ./a.out
0.860 s
 $ ./sub pians.txt output.txt
-0.39999999996816507571e-99999

結果

• 2万桁は 2.8秒 → 0.03秒
• 5万桁は21.23秒 → 0.22秒
• 10万桁は0.860秒
• ... 目標値は 100万桁で10秒

バリ糞速くなりましたねwww 100倍速くなったwww

100万桁, そして1000万桁の世界へ

このまま100万桁, そして目標の1000万桁まで行きましょう. pians.txtを100万桁まで書き*17, sub.cのdigitsを1000000にします.

20万桁 (pi.cのdigitsを200000に)

 $ gcc pi.c -lgmp && ./a.out
3.650 s
 $ ./sub pians.txt output.txt
-0.14797927213959375309e-199999

40万桁

 $ gcc pi.c -lgmp && ./a.out
17.820 s
 $ ./sub pians.txt output.txt
-0.31728225152177753651e-399999

うー... 10秒超えました...

どうやらこのままでは100万桁はまだまだ遠いようです. プログラムの改善が必要です.

ここでBellardさんの記事を見返します. 実はこの記事を最初から手に持っていたら分かるんですが, ｢Binary Splitting method｣を使いなさい, という事なんですね. え, Bina...何それ? 今から書きますよ.

もう一度, 数式に戻ります. さっきはこの漸化式を使いました.

しかし, これはだいたい, (桁の大きな数) * (桁の小さな数) の計算ばかりになるので, 遅いのです. これを次のようにします.

について

さっき使った漸化式は, これのと全く同じです. しかし, Binary Splitting methodではとします. 掛け算の数はおそらく変わりません. しかし, 桁の数がめちゃくちゃ違うようなアンバランスな掛け算が減るので, 効率が良くなるのです.

Binary Splitting method によるChudnovsky実装

再帰でガーーーーッって行く感じ, 割と直感的に書きました. 以下のコードがこのエントリーの最終的なコードです. ベンチを取るときはこのコードを使用してください.

/* Calculate pi based on Chudnovsky algorithm (& Binary Splitting method), using GMP */
/* [1] Computation of 2700 billion decimal digits of Pi using a Desktop Computer,
       Fabrice Bellard, Feb 11 2010 (4th revision),
       http://bellard.org/pi/pi2700e9/pipcrecord.pdf */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>

mpz_t A, B, C, C3over24;

void computePQT(mpz_t P, mpz_t Q, mpz_t T, int n1, int n2) {
  int m;
  if (n1 + 1 == n2) {
    mpz_t a;
    mpz_init(a);
    /* P(n2-1, n2) = - (2 * n2 - 1) * (6 * n2 - 5) * (6 * n2 - 1) */
    mpz_set_ui(P, (2 * n2 - 1));
    mpz_mul_ui(P, P, (6 * n2 - 5));
    mpz_mul_ui(P, P, (6 * n2 - 1));
    mpz_neg(P, P);
    /* Q(n2-1, n2) = C^3 * n2^3 / 24 */
    mpz_mul_ui(Q, C3over24, n2);
    mpz_mul_ui(Q, Q, n2);
    mpz_mul_ui(Q, Q, n2);
    /* a_n2 = (A + B * n2) */
    mpz_mul_ui(a, B, n2);
    mpz_add(a, a, A);
    /* T(n2-1, n2) = a_n2 * P(n2-1, n2) */
    mpz_mul(T, a, P);
    mpz_clear(a);
  } else {
    mpz_t P1, Q1, T1, P2, Q2, T2;
    mpz_init(P1);
    mpz_init(Q1);
    mpz_init(T1);
    mpz_init(P2);
    mpz_init(Q2);
    mpz_init(T2);
    m = (n1 + n2) / 2;
    computePQT(P1, Q1, T1, n1, m);
    computePQT(P2, Q2, T2, m, n2);
    /* P = P1 * P2 */
    mpz_mul(P, P1, P2);
    /* Q = Q1 * Q2 */
    mpz_mul(Q, Q1, Q2);
    /* T = T1 * Q2 + P1 * T2 */
    mpz_mul(T, T1, Q2);
    mpz_mul(P1, P1, T2);
    mpz_add(T, T, P1);
    mpz_clear(P1);
    mpz_clear(Q1);
    mpz_clear(T1);
    mpz_clear(P2);
    mpz_clear(Q2);
    mpz_clear(T2);
  }
}

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

  clock_t start, end;
  FILE* fp;

  int digits = 10000;
  int prec = digits * log2(10);
  int n = digits / 14;

  start = clock();

  mpf_t pi, temp;
  mpz_t P, Q, T;

  /* Initialize GMP numbers */
  mpf_set_default_prec(prec);
  mpf_init(pi);
  mpf_init(temp);
  mpz_init(P);
  mpz_init(Q);
  mpz_init(T);
  mpz_init(A);
  mpz_init(B);
  mpz_init(C);
  mpz_init(C3over24);

  /* Assignment */
  mpz_set_str(A, "13591409", 10);
  mpz_set_str(B, "545140134", 10);
  mpz_set_str(C, "640320", 10);
  mpz_mul(C3over24, C, C);
  mpz_mul(C3over24, C3over24, C);
  mpz_div_ui(C3over24, C3over24, 24);

  computePQT(P, Q, T, 0, n);

  /* pi = C ^ (1 / 2)     */
  mpf_set_z(temp, C);
  mpf_sqrt(pi, temp);
  /*      * C             */
  mpf_mul(pi, pi, temp);
  /*      * Q             */
  mpf_set_z(temp, Q);
  mpf_mul(pi, pi, temp);
  /*      / (T + A * Q)   */
  mpz_mul(Q, A, Q);
  mpz_add(Q, Q, T);
  mpf_set_z(temp, Q);
  mpf_div(pi, pi, temp);
  /*      / 12            */
  mpf_div_ui(pi, pi, 12);

  /* mpf_out_str(stdout, 10, digits, pi); */
  fp = fopen("output.txt", "w");
  if (fp == NULL) {
    printf("couldn't open output.txt");
    exit(EXIT_FAILURE);
  }
  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);
  mpf_out_str(fp, 10, digits, pi);

  mpf_clear(pi);
  mpf_clear(temp);
  mpz_clear(P);
  mpz_clear(Q);
  mpz_clear(T);
  mpz_clear(A);
  mpz_clear(B);
  mpz_clear(C);
  mpz_clear(C3over24);

  end = clock();
  printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

  return 0;
}

で, 実行結果

 $ # 20万桁
 $ gcc pi.c -lgmp && ./a.out
0.320 s
0.400 s
 $ ./sub pians.txt output.txt
-0.14797927213959375309e-199999
 $ # 40万桁
 $ gcc pi.c -lgmp && ./a.out
0.780 s
1.000 s
 $ ./sub pians.txt output.txt
-0.31728225152177753651e-399999
 $ # 100万桁
 $ gcc pi.c -lgmp && ./a.out
2.510 s
3.220 s
 $ ./sub pians.txt output.txt
0.13092756283208453158e-999999

結果

• 20万桁は 3.65秒 → 0.40秒 (10倍速くなったw)
• 40万桁は 17.82秒 → 1.00秒 (これも10倍www)
• 60万桁は 1.65秒
• 80万桁は 2.41秒
• 100万桁は 3.22秒
• → →目標値の 100万桁で10秒 を余裕をもってクリア!!!

(∩´∀｀)∩ﾜｰｲ 円周率100万桁を3.22秒で計算できたよ! ヽ(´ー｀)ノﾊﾞﾝｻﾞｰｲ

1000万桁, そして1億桁を目指して

1000万桁まで行きましょう. 今回は行くことが出来ました.

• 200万桁 7.81秒 差: -0.38782687488202153327e-2000000
• 400万桁 18.89秒 差: 0.40845608946106963042e-4000000
• 800万桁 44.82秒 差: 0.46321418078726345026e-7999999
• 1000万桁 59.27秒 差: -0.27408486638735177309e-9999999

1000万桁で一分程度(計算は47.33秒です. 59.27秒っていうのはファイルへの書き込みを含めてのタイムです). 当初の目標である1000万桁は達成されました!!!

人間ここまで来ると欲望が溢れでてくるわけです. 次は1億桁です. ここで, 最適化を試みました.

• computePQTのなかで, P1, ... T2と六つも変数を作っては破棄してをしていますが, P2, Q1, T1をそれぞれ引数のP, Q, Tで代用できます. しかしこれをやってもあまり計算は速くならなかったので(変数名が分けがわからなくなるし)辞めました.
• 条件分岐で a = (-1)^ ... の符号を反転させてるけど, 実はaではなくてpの符号を｢常に｣反転させても同じ結果になります. これ, ちょっとだけ速くなりました. (上のコードはこの変更を反映させています)
• を計算する所で, を後で掛けたり先に掛けたりしましたが, 変わりませんでした.
• Bellardさんの記事2.2.1によると, とは, これらの最大公約数で割ってもOKです. ところが, これをmpz_gcd, mpz_divで実装すると, 逆に遅くなりました. modの計算は速いとしても, 割り算2回をいちいちするのが遅くなる原因かもしれません.

　

• Bellardさんの記事2.2.2によると, , を見直すことで更に高速化できると書いてあります. ところが, これを実装してもあんまり速くなりませんでした.

結局3時間くらい高速化しようと頑張っていましたが, Bellardさんの記事に書いてある方法も, あまり速くはなりませんでした. 多分, 彼らのやっているような｢何兆桁｣のレベルになると, 割り算をすることで増えるコストよりも, 最大公約数で割って桁数を下げることで減るコストの方が支配的になって, 速くなるんだと思います.

結局高速化は諦めて, 先ほどのコードを最終形としてどんどん計算しました. で, 最後には1億桁達成へ!!!

桁数	計算時間	トータル時間	実際の値との差
10,000	0.00s	0.00s	-0.14332772033801142172e-9999
100,000	0.12s	0.16s	-0.35873997562031545622e-99999
200,000	0.32s	0.41s	-0.14797927213959375309e-199999
400,000	0.78s	1.00s	-0.31728225152177753651e-399999
800,000	1.90s	2.42s	-0.2680842256059386162e-799999
1,000,000	2.65s	3.38s	0.13092756283208453158e-999999
2,000,000	6.12s	7.77s	-0.38782687488202153327e-2000000
4,000,000	14.89s	18.79s	0.40845608946106963042e-4000000
8,000,000	35.76s	44.81s	0.46321418078726345026e-7999999
10,000,000	46.83s	58.86s	-0.27408486638735177309e-9999999
16,777,216	92.05s	114.05s	0.18253909916599945649e-16777215
20,000,000	110.65s	138.16s	0.41731393365811761333e-19999999
40,000,000	264.63s	328.26s	-0.208525002182786562e-39999999
60,000,000	429.02s	533.13s	0.13663586443522920956e-59999999
80,000,000	621.41s	765.73s	0.28594903035388918799e-79999999
100,000,000	807.20s	995.11s	0.2e-99999999 (ただし, pians.txtが1億桁)

ヨッシャー!!! ∩(＞◡＜*)∩ 一億桁達成や!!!

最後の記録は, pians.txtの方が1億桁までしか書いていなかったので, ちょっと残念... ですが, 1億桁行ったって言っていいよね...!!!

トータル時間 = 計算時間 + ファイル出力時間 です. 検証の時間は含みません. ファイル出力にえらい時間かかっているなぁという印象を持たれるかもしれませんが, この過程では, 数字の16進数表現から10進数表現への変換もありますので, こんなもんでしょう.

大事なのはオーダーですね. まずは時間をグラフにしてみました. 以下では, ファイル出力を除いた, 計算時間のみを扱っています. まずは生の時間.

緩やかに増えていきます. 二次関数ではなさそうですが, 一次関数よりはでかそうです. 実のところ, ChudnovskyアルゴリズムはO(n log(n)^3)らしいですが*18, どうでしょうか? 計算時間をn log(n)^3で割ってみます.

うーん... n log(n)^4で割ってみますと,

これかなぁ... O(n log(n)^4)っぽいです. なんでこうなったかは分かりません.

この次のステップへ

なんとか一億桁まで割と現実的な時間で計算をし, 検証まで行いましたが, 以下の点がダメダメです.

• メモリーのこととかスワップのこととかを一切考えていない. コンパイラとOS任せ. これ以上, 桁を伸ばそうとすると, 数字がメモリーに乗らなくなるので, おそらくバイナリーデータのままファイルに一時保存したり再度読みだしたりしながら計算するんだと思う. よく分からない.
• 検証方法がダメ. 今回は, 1億桁の円周率をどこからか*19取ってきて, これを一行のテキストファイルに直して, 自分のプログラムの出力をファイルに書きこんで, 更に2つのファイルを読み込んで, 引き算しました. ですが, 実際の円周率計算では, 未知の桁数を計算します. よって, こんなのはチートです. でもよく分からない... (そして, ./sub の実行時間を一切書いていませんが, 実は桁数が増えると割と時間がかかります. ですから目標桁数+100桁くらいのファイルを作っておいて, それと比較しています.)
• もうちょっと速くならないかなぁ... gcdを取るという方法もアリですし, まだ何かありそうです. 重要なのは, コードを見つめるのではなく, 数式を見つめることです. GMPはとても良くできているし, テストも多く書かれているでしょうから, 掛け算を細々と最適化した所で殆ど変わりません. 数式を見つめ, 数字を小さくすること, これに尽きますね.
• というかy-cruncher速すぎワロタwww このソフトのサイトのベンチマーク, 自分の結果と一桁違いまっせwww Windows用のソフトだから, 自分のPCで確かめられないのが残念だわwww

結論

• Chudnovskyアルゴリズムは素晴らしいです. 他の方法は知らないけど, Chudnovskyの速さは皆さん(円周率計算トップレベルの方々)のお墨付きですからね.
• 私は円周率計算の素人だなぁと思いました. まだまだ技術的に難しい点があります.
• 滅茶苦茶楽しかった. アルゴリズムを書き換えたらどんどん速くなっていくのが超楽しかった. それのみならず, 1万桁まで上手く動くコードが書けたら, そのコードで10万桁でも100万桁でも1000万桁でもきちんと動くっていうのが楽しかった. 楽しいことは一番だ!
• ちょっとこの後を頑張るのは辞めようかなぁと. 近藤さんのUstreamを見たから分かるけど, 彼らのやっていることは｢10兆桁｣です. 私が計算した｢1億桁｣の10万倍ですか... とてつもないです. HDDとの戦いになります. それに, 3年前に購入したへっぽこPCのCPUで, しかもシングルコアで戦っている私とはレベルが違います. (96GBのメモリーだとか, 2TBのHDDが20個だとか... *20 ) 桁が違いすぎて鼻血が出そうですね. プログラミングのレベルも相当なものだと思います. 尊敬します.

みなさんも円周率の計算, やってみてはいかがでしょうか. そして私を超えて行って下さい!!!

追記: 他のプログラムとの計算速度の比較 [2012/03/05]

せっかく円周率への熱がまだ残っているうちに, 他の人が書いたプログラムとの比較をしておきます.

番号	1	2	3	4
プログラム名	pi.c	QuickPi v4.5	PI Calculation Program	Hironobu SUZUKIさんのプログラム
アルゴリズム	Chudnovsky	Chudnovskys他	Gauss AGM	Gauss AGM
実装した人	itchyny	Steve Pagliarulo	Takuya OOURA	Hironobu SUZUKI
URL	http://d.hatena.ne.jp/itchyny/20120304/1330870932	http://members.shaw.ca/francislyster/pi/pi.html	http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html	http://h2np.net/pi/
備考	GMP/Bellardさんの記事に基づいてプログラムを書いた, 時間はファイル出力含む	他のアルゴリズムも使用できる. Dual coreを思いっきり使う	pi_ca使用/2の冪乗っぽい数字の桁数しか計算できない	GMP/ソースコードはurl参照, 1000万桁の計算でアウト?
10,000桁	0.00s	0.00s	0s	0.020s
	(-0.143e-9999)	(-0.566e-10000)	(-0.790e-12289)	(-0.38e-10009)
100,000桁	0.16s	0.07s	0s	0.47s
	(-0.358e-99999)	(-0.413e-100000)	(0.528e-163840)	(-0.203e-100009)
1,000,000桁	3.38s	1.20s	5s	7.56s
	(0.131e-999999)	(-0.309e-1000000)	(-0.556e-1310720)	(0.321e-1000009)
10,000,000桁(1000万)	58.86s	19.42s	72s	116.07s
	(-0.274e-9999999)	(-0.259e-10000000)	(-0.291e-8388603)	(0.63e-1430651)
16,777,216桁	114.05s	36.45s	161s
	(0.183e-16777215)	(-0.825e-16777216)	(-0.122e-16777211)
100,000,000桁(1億)	995.11s	346.69s	3445s (len of FFT=20000000)
	(0.2e-99999999)	(0e-100000000)	(1e-134217727)

番号	5	6	7	8
プログラム名	pyopyopyoさんのプログラム	ramanujan.c	ghc-chudnovsky.c	PI.hs
アルゴリズム	Square AGM	Borwein's algorithm	Chudnovsky	Chudnovsky
実装した人	pyopyopyo	itchyny	Hanhong Xue	tanakh
URL	http://d.hatena.ne.jp/pyopyopyo/20090308/p1	http://mailsrv.nara-edu.ac.jp/~asait/c_program/sample0/pi_source.htm	http://gmplib.org/pi-with-gmp.html	https://github.com/tanakh/pi, http://tanakh.jp/posts/2012-03-08-pi.html
備考	GMP (gmpxx.h)/ソースコードはurl参照	GMP使用	GMP公式	Haskellによる
10,000桁	0.020s	0.02s (loop:5)	0.004s	0.02s
	(-0.143e-9999)	(-0.143e-9999)	(-0.333e-10001)
100,000桁	0.74s	0.70s (loop:8)	0.108s	0.35s
	(-0.359e-99999)	(-0.631e-99994)	(0.126e-100001)
1,000,000桁	13.50s	13.45s(loop:10)	1.916s	6.29s
	(0.131e-999999)	(0.131e-999999)	(0.928e-1000002)
10,000,000桁(1000万)	239.41s	247.40s (loop:13)	34.706s	121s
	(-0.274e-9999999)	(-0.134e-9999996)	(-0.408e-10000001)
16,777,216桁	433.19s		64.584s	220s
	(0.183e-16777215)		(0.254e-16777217)
100,000,000桁(1億)		3840.52s (loop:14)	575.039s	507s
		(-0.28e99999998)	1e-100000000

ざっと横に見てみますと, QuickPiの速さがピカイチです. 自分のコードはそれに比べて3倍遅いですが, オーダーは似たようなものです. AGMによる実装は軒並遅いなぁっていう感じがしますね. 一番右のアルゴリズムについては後で言及します. それぞれのプログラムについて軽くコメントしておきます.

1. このエントリーで示したプログラムそのものです. 1000万桁は一分程度, 1億桁は16分程度でした.
2. Steve Pagliaruloさんが配布されていらっしゃるプログラム. これはソースコードはありません. Chudnovskyで, 1000万桁20秒, 1億桁が6分弱と, 私のプログラムのおよそ1/3のタイムで実行できました.
3. Takuya OOURAさんが配布されていらっしゃるプログラム. ソースコードで配布されており, makeでコンパイルするようになっています. 詳細は見ていませんが, FFTのパッケージのようです. 円周率のアルゴリズムはAGM (Arithmetic Geometric Mean)です. 実行時間は1000万桁で一分を超えています.
4. Hironobu SUZUKIさんのホームページにあったコードです. 少し書き換えて実行しています. これもAGMで, 1000万桁で二分程度かかっています.
5. pyopyopyoさんのブログのソースコードを若干改変して使用させて頂きました. アルゴリズムはAGM(前2つと何が違うか知りません)で, 1000万桁の計算で4分弱かかっています. [追記:date=2012/03/09]前2つと同じアルゴリズムであるとのことです. id:periaさん, ありがとうございました.[/追記]
6. Ramanujanによるアルゴリズムらしいです. 実装は私が行いました. ソースコードは後ろに載せています. 参考URLの論文は後日読みます. アルゴリズムの名称は知りません. 参考URLを教えてくださった @kagakuma さんに感謝します. [追記:date=2012/03/09]Ramanujanは関係なくて, Borwein's algorithm(Borweinの4 次の公式)と言うそうです*21. 1ステップ進むごとに4倍の精度になるんですね. id:periaさん, ありがとうございました.[/追記]
   

   15回の計算で円周率を10億桁以上求める方法 URL 小学生の時からずっと気になってるんだけど、誰か試してみてくれませんか？プログラミングとか駆使したらいけるのかな。

   2012-03-04 22:24:40 via Tween

7. GMPのサイトで円周率のベンチマークとして掲載されているプログラム*22. やはりChudnovskyアルゴリズムで, 独自の高速化の工夫が光りました.
8. id:tanakhさんによるHaskellでの実装. やはりChudnovskyで, これもなかなか速いです. ただ, 1億桁の結果には少し疑問が残ります. リダイレクトでoutput.txtに書きこんだのですが, 書き込みがうまく行ってなかったからです. 1000万桁の時はこのような事はありませんでした. 原因不明...

上記以外の, 円周率の計算するプログラムとして有名なものは以下のものがありますが, 何れも実行できませんでした. Windowsをお使いの人は, y-cruncherはやってみて価値があると思います.

• TachusPI: 手持ちのPCが32bitのため, 実行できず. ソースコード配布なし. アルゴリズム: Chudnovsky, 作者: Fabrice Bellard, URL: http://bellard.org/pi/pi2700e9/tpi.html
• y-cruncher: Windows向けのため, 実行できず. ソースコード配布なし. アルゴリズム: Chudnovsky, 作者: Alexander J. Yee URL: http://www.numberworld.org/y-cruncher/

Hironobu SUZUKIさんのプログラムは若干変更して実行しています.

7a8
> #include <stdlib.h>
15a17,18
>   clock_t start, end;
>   FILE* fp;
17a21
>   mpz_t temp, temp2;
19c23,24
<   dprec = 1000000L;		/* decimal precision */
---
>   start = clock();
>   dprec = 10000000L;		/* decimal precision */
34a40,41
>   mpz_init (temp);
>   mpz_init (temp2);
41a49
>   mpz_set_ui (temp2, 2);
53c61,62
<       mpf_pow_ui (t2, c2, k);	/* 2^k */
---
>       mpz_pow_ui (temp, temp2, k);	/* 2^k */
>       mpf_set_z (t2, temp);
59c68
<   mpf_pow_ui (t2, t1, 2);
---
>   mpf_mul (t2, t1, t1);
62c71,78
<   mpf_out_str (stdout, 10, dprec + 10, A);
---
>   fp = fopen("output.txt", "w");
>   if (fp == NULL) {
>     printf("couldn't open output.txt");
>     exit(EXIT_FAILURE);
>   }
>   mpf_out_str(fp, 10, dprec + 10, A);
>   end = clock();
>   printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

pyopyopyoさんのプログラムも変更させていただきました.

5a6,7
> #include <stdio.h>
> #include <stdlib.h>
11c13
< #define  N  100*100*100
---
> #define  N  10000000
17a20,22
>      clock_t start, end;
>      FILE* fp;
>      start = clock();
54a60
>      }
60,70c66,69
<      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 <<" ";
---
>      fp = fopen("output.txt", "w");
>      if (fp == NULL) {
>        printf("couldn't open output.txt");
>        exit(EXIT_FAILURE);
71a71,73
>      mpf_out_str(fp, 10, N, tmp.get_mpf_t());
>      end = clock();
>      printf("%.3f s\n",(double)(end - start) / CLOCKS_PER_SEC);

一番右のRamanujanによるアルゴリズムについてここで少し言及しておきます*23. ただし, 私は論文はまだ確認していません. 数列とを次のように定義します. [追記:date=2012/03/09]Ramanujanは関係ありません. Borwein's algorithmと言うそうです.[/追記]

この初項と漸化式に基づいて計算していくと, 大きなで

となるんですって. ほんまかいなと思う前にキーボードを叩きコードを書いて実装してみたら, 本当にそうなったからほんまなんでしょう.

やはりGMPを用いた実装です.

/* Calculate pi using Borwein's algorithm */
/* [1] J.Borwein and P.Borwein,  Ramanujan and Pi,
       IEEE, Science and Applications, Supercomputing 88, Volume II, 117-128,
       1988 */
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
int main(int argc, char *argv[]) {

  clock_t start, end;
  FILE* fp;
  int digits = 100000;
  int prec = digits * log2(10);
  int loop = 8;
  int i;

  start = clock();
  mpf_t y, a, temp, twopower;
  mpf_set_default_prec(prec);
  mpf_init(y);
  mpf_init(a);
  mpf_init(temp);
  mpf_init(twopower);

  /* y0 = 2 ^ (1 / 2) - 1; a0 = 6 - 4 * 2 ^ (1 / 2); */
  mpf_set_ui(y, 2);
  mpf_sqrt(y, y);
  mpf_set(a, y);
  mpf_sub_ui(y, y, 1);
  mpf_mul_ui(a, a, 4);
  mpf_sub_ui(a, a, 6);
  mpf_neg(a, a);
  mpf_set_ui(twopower, 2);

  for (i = 1; i <= loop; i++) {

    /* y[i] = (1 - (1 - y[i-1] ^ 4) ^ (1 / 4)) / (1 + (1 - y[i-1] ^ 4) ^ (1 / 4)) */
    mpf_mul(y, y, y);
    mpf_mul(y, y, y);
    mpf_sub_ui(y, y, 1);
    mpf_neg(y, y);
    mpf_sqrt(y, y);
    mpf_sqrt(y, y);
    mpf_add_ui(y, y, 1);
    mpf_ui_div(y, 2, y);
    mpf_sub_ui(y, y, 1);

    /* a[i] = (1 + y[i]) ^ 4 a[i-1] - 2 ^ (2 * i + 1) * y[i] * (1 + y[i] + y[i] ^ 2) */
    mpf_add_ui(temp, y, 1);
    mpf_mul(temp, temp, temp);
    mpf_mul(temp, temp, temp);
    mpf_mul(a, a, temp);
    mpf_mul_ui(twopower, twopower, 4);
    mpf_mul(temp, y, y);
    mpf_add(temp, temp, y);
    mpf_add_ui(temp, temp, 1);
    mpf_mul(temp, temp, y);
    mpf_mul(temp, temp, twopower);
    mpf_sub(a, a, temp);

  }

  /* pi = 1 / a */
  mpf_ui_div(a, 1, a);

  fp = fopen("output.txt", "w");
  if (fp == NULL) {
   printf("couldn't open output.txt");
   exit(EXIT_FAILURE);
  }
  mpf_out_str(fp, 10, digits, a);
  end = clock();
  printf("%.3lf s\n",(double)(end - start) / CLOCKS_PER_SEC);

  return 0;
}

loopとdigitsを適当に変えながら出力を見てください. このアルゴリズム, ホントワケが分かりませんが, きちんと動くことは動きます. ただ, さっきの表を見てください. 1000万桁で4分程度なんですが, 1億桁を計算したら3840秒, これって64分, 1時間超えてるんです. 14回ループして1億桁たぶん超えてますけどね.

言いたいことは, ｢ループの回数が少ないから速いというわけではない｣という事です. なぜでしょうか. もう言うまでもないことかもしれませんね. 浮動小数点のせいです. 1億桁の計算をしようと思ったら, 浮動小数点の変数ははじめから1億桁です(バイナリーはこれのlog2(10)倍). それらを1億桁の精度を持って掛け算し, 割り算し, しかも平方根を取ったりしています. pyopyopyoさんのプログラムもHironobu SUZUKIさんのプログラムもそうです. ループの中が浮動小数点演算という時点で, 計算速度が良くなるはずがないのです. 割り算のコストもたぶんやばいことになっています.

それに対して, Chudnovskyが圧倒的に速いのは, このアルゴリズムはループの中が全て整数の演算(しかも掛け算か足し算, 符号反転のみ, 割り算はありません)に収まっているからです. 最初に浮動小数点で書いておいて, それを整数に書きなおしたら｢バリ糞速く｣なりましたよね.しかも, それは, このエントリーで試行錯誤したように, "Binary Splitting method"によって更に劇的に速くなる. この方法によって整数たちはお互いにどんどん小さくなり, 同等の長さ同士の掛け算になるんです.

RamanujanのBorweinのアルゴリズムの参考URLではi=15に達したら20億桁まで計算できると書いてあります. しかし, それはこのアルゴリズムが良いという事では決してないことを, ここで強く述べたいと思います. 漸化式見たら, ループ毎に平方根2回としかも割り算付きですよ. それを20億桁の数字でやろうと思ったら... (もういいですね... 一時間もかかったからちょっと腹がたってしまって...) あと, 一時間もしたらclock_tが桁溢れします. 出力されるマイナスの数から逆演算しました. というか, そもそもprintfのフォーマット指定が間違ってるような...? まぁええか.

追記: 一億桁を｢きちんと｣超えておく [2012/03/05]

一億桁, 超えてないんじゃない!って突っ込まれる気がしたので, きちんと超えておきました. 計算時間は13分36秒, トータル時間は16分45秒です.

 $ # int digits = 1000100; で計算しました.
 $ gcc pi.c -lgmp && ./a.out
816.210 s
1005.390 s
 $ hexdump -C output.txt | head -n 10
00000000  30 2e 33 31 34 31 35 39  32 36 35 33 35 38 39 37  |0.31415926535897|
00000010  39 33 32 33 38 34 36 32  36 34 33 33 38 33 32 37  |9323846264338327|
00000020  39 35 30 32 38 38 34 31  39 37 31 36 39 33 39 39  |9502884197169399|
00000030  33 37 35 31 30 35 38 32  30 39 37 34 39 34 34 35  |3751058209749445|
00000040  39 32 33 30 37 38 31 36  34 30 36 32 38 36 32 30  |9230781640628620|
00000050  38 39 39 38 36 32 38 30  33 34 38 32 35 33 34 32  |8998628034825342|
00000060  31 31 37 30 36 37 39 38  32 31 34 38 30 38 36 35  |1170679821480865|
00000070  31 33 32 38 32 33 30 36  36 34 37 30 39 33 38 34  |1328230664709384|
00000080  34 36 30 39 35 35 30 35  38 32 32 33 31 37 32 35  |4609550582231725|
00000090  33 35 39 34 30 38 31 32  38 34 38 31 31 31 37 34  |3594081284811174|
 $ hexdump -C output.txt | tail -n 10
05f5e0e0  39 38 31 36 31 31 36 38  33 31 33 39 33 37 35 31  |9816116831393751|
05f5e0f0  34 39 37 30 35 38 31 31  32 30 31 38 37 37 35 31  |4970581120187751|
05f5e100  35 39 32 32 31 35 30 35  38 38 30 39 35 37 38 33  |5922150588095783|
05f5e110  32 37 39 36 33 34 38 37  33 30 39 35 31 33 35 32  |2796348730951352|
05f5e120  38 34 39 31 31 30 33 33  34 31 37 39 37 35 37 32  |8491103341797572|
05f5e130  30 31 32 35 38 38 33 34  30 36 32 31 33 36 39 30  |0125883406213690|
05f5e140  35 34 32 32 39 35 38 33  38 37 38 39 34 36 30 37  |5422958387894607|
05f5e150  31 34 32 34 38 35 35 39  37 32 32 31 30 30 38 34  |1424855972210084|
05f5e160  38 31 35 36 36 31 65 31                           |815661e1|
05f5e168

0x05f5e100がぴったり1億桁です. つまり, このファイルの1億1番目(1足したよ!)から 5 9 2 2 1 5 ...となっています. しかし最初の"0.3"があるので, 三文字ずらして, 円周率小数点以下1億1番目の数は2, 次は1, そして5... となります. 10文字ずつにすると,

2150588095 7832796348 7309513528 4911033417 9757201258 8340621369 0542295838 7894607142 4855972210 084815661e1

これと, One billion digits*24に100000001をsubmitした結果を突き合わせます.

2150588095 7832796348 7309513528 4911033417 9757201258 8340621369 0542295838 7894607142 4855972210 0848156605 8666322078 8245761718

1億1桁 1億11桁, 21, ... 1億91 (0), 92 (8), 93 (4) 94 (8) 95 (1) 96 (5) 97 (6) 98 (6) 99 (1 vs 0)

という訳で, 1億98桁まで正確に求められたことになります. これで安心して｢1億超えた!!!｣って言えますね. 安心安心.

追記: GMPのC++インターフェースを使ってみる [2012/03/05]

GMPとC++で更にコードが書きやすくなります. (C++は普段全然書かないので, おかしなとこがあるかもです.)

/* Calculate pi based on Chudnovsky algorithm (& Binary Splitting method), using GMP */
/* [1] Computation of 2700 billion decimal digits of Pi using a Desktop Computer,
       Fabrice Bellard, Feb 11 2010 (4th revision),
       http://bellard.org/pi/pi2700e9/pipcrecord.pdf */
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <fstream>
#include <iostream>
#include <gmpxx.h>

mpz_class A, B, C, C3over24;
class PQT {
  public:
    mpz_class P;
    mpz_class Q;
    mpz_class T;
};

PQT computePQT(int n1, int n2) {
  int m;
  PQT ans;
  if (n1 + 1 == n2) {
    ans.P =  - (2 * n2 - 1);
    ans.P *= (6 * n2 - 5);
    ans.P *= (6 * n2 - 1);
    ans.Q = C3over24 * n2 * n2 * n2;
    ans.T = (A + B * n2) * ans.P;
  } else {
    m = (n1 + n2) / 2;
    PQT PQT1 = computePQT(n1, m);
    PQT PQT2 = computePQT(m, n2);
    ans.P = PQT1.P * PQT2.P;
    ans.Q = PQT1.Q * PQT2.Q;
    ans.T = PQT1.T * PQT2.Q + PQT1.P * PQT2.T;
  }
  return ans;
}

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

  clock_t start, end;
  int digits = 1000000;
  int prec = digits * log2(10);
  int n = digits / 14;
  start = clock();
  mpf_class pi (0, prec);
  A = "13591409";
  B = "545140134";
  C = "640320";
  C3over24 = C * C * C / 24;
  PQT PQT = computePQT(0, n);
  pi = C;
  pi = sqrt(pi);
  pi *= C * PQT.Q;
  pi /= (PQT.T + A * PQT.Q) * 12;

  end = clock();
  std::cout << (double)((end - start) / CLOCKS_PER_SEC) << "s" << std::endl;

  std::ofstream ofs ("output.txt");
  ofs.precision(digits);
  ofs << pi << std::endl;
  end = clock();
  std::cout << (double)((end - start) / CLOCKS_PER_SEC) << "s" << std::endl;

  return 0;
}

 $ g++ pi.cpp -lgmpxx -lgmp && ./a.out
2s
3s

 $ ./sub output.txt pians.txt
0.413e-100000

100万桁で3秒, 1000万桁で48秒(出力あわせて60秒)でした. C言語の場合と速度はあまり変わりません.

一つだけ注意があります.

  pi = C;
  pi = sqrt(pi);
  pi *= C * PQT.Q;
  pi /= (PQT.T + A * PQT.Q) * 12;

ここですが,

  pi = sqrt(C) * C * PQT.Q / (PQT.T + A * PQT.Q) / 12;

と書くとうまく行きません. 理由はもう言わずもがなですよね. C++のインターフェースだったら, ここら辺ももっともっと書きやすくなてもいいのになぁと思いました. 書き易さと落とし穴で, なんか中途半端な感じがして私は嫌いだなぁ.

追記: 他の方の実行結果[2012/03/05]

このエントリーに載せたプログラムの実行結果(円周率ではなくて, 出力含めた計算時間)を報告していただければ, ここに載せることにします. Twitterでもブログでもブクマコメでも, お待ちしております.

• ファイル出力を含めた時間をお願いします.
• CPU, OSもあわせてお願いします.
• プログラムの実行に伴う如何なる破損についても責任を負いかねます.
• 私のパソコンで1億が16分ですので, その程度の計算時間は覚悟してください.
• 1億桁でのoutput.txtは95MBとなります. 普通のテキストエディタでは開かないようにしてください.

お名前	PC	1000万桁	一億桁	十億桁
鏡 弘道さん @kagami_hr	2.7GHz MacBook Pro	16.5s	261s
		@itchyny とりあえず一千万桁でやって itchyny さんのプログラムが出力込みで 16.5 秒. CLN が 14.8 + 7.4 秒位です。今一億桁やっています。 2012-03-05 15:25:46 via YoruFukurou to @itchyny	@itchyny 実時間で 4m21.114s = 261 秒。速いです。 2012-03-05 15:32:30 via YoruFukurou to @itchyny
いちにぃ @itchyny	Intel(R) Core(TM)2 Duo CPU P8600 @2.40GHz Ubuntu	58.86s	995.11s
			やったぁ！ 円周率一億桁いったぁ！ (計算時間: 807秒) 2012-03-04 22:23:27 via ぼくは二次元の女の子が好きだ (※注意:これは出力の時間除いてのタイムですね)
id:nitobeさん	3.4G Intel Core i7-2600, debian64, VirtualBox, Windows 7		227.850s
id:nitobeさん	3.4G Intel Core i7-2600, debian64, VirtualBox, Windows 7, gcc-4.4.5, gmp-5.1.0		200.16s	3059.55s

ツイートする

Permalink | コメント(8) | トラックバック(2) | 23:22