[NEW] 最新のgmp-4.3.0を利用し、AMD64 / Intel64とも、さらに高速化しました。[20090421]

GNU/Linux上で円周率の計算をおこなう

はじめに

そもそもは、円周率計算のプログラムをケースにしてGMPの能力の限界はどれくらいかということを試していたのだが、手段が目的化してしまい、そのなれの果てがこの結果である。

まずは小数点以下100万桁までを求めてみる

準備

GMP (GNU Multi-Precision)ライブラリを使うが、ディストリビューションのデフォルトでついてくるライブラリは遅い。正確に言うと、使っているハードウェアに最適化していない汎用の環境を前提にライブラリはコンパイルされているため、さらに最適化することができる。

GMPのソースコードを入手する ( http://gmplib.org/ )
ソースコード展開後、コンフィグレーション、コンパイル、インストールを行う

デフォルトは/usr/local/以下にインストールされる。筆者のAthlon 1GHzを使っている環境ではDebian のデフォルトのgmpと、サイトで独自にコンパイルしたものでは約25%程度の速度の違いがあった。

GMPに関しては、さらにチューンアップができるが、筆者が試行錯誤してアップしたのは2-3パーセントであった。苦労した割にはあまり成果はあがらなかった。当然ながら、これは使用している環境に強く依存しているので、ある環境上では、十分な成果があるかも知れない。試してみたい人はgmp/tuneの下にチューンアップのためのツールがあるので、それを個々の環境で試して欲しい。

円周率を求めるプログラムを書く

円周率を求める計算方法は数々あるが、高速に収束するGauss AGM法を用いる。下記のプログラムにはSchonhage方式のGauss AGM法を使っている。サンプルのコードは下記のようなものである。ファイルはpi.cである。

Gauss AGM法は絶対誤差 (Absolute Error)が残るが、プログラム中では収束させるためのループが十分なので、出力中の桁の中には誤差が現われない。計算には100万桁にプラス数桁の余裕を持たせて計算させているのは、別の検証プログラムと誤差を比較する時に使っていた名残である。

下記の記述はgmp-3.0の頃の記述である。gmp-4.0になってからは、出力が高速化されているので、この記述はあてはまらない。(2006年8月17日)

100桁分の計算は高速であるが、出力には時間がかかる。このプログラムでもっとも時間がかかるのは、最後の行のmpf_out_str ()関数である。筆者が使っているAthlon 1G, Pentium 4 1.6G のどちらのCPU上でも、100万桁の計算には50 秒未満だが、桁出力には100秒以上かかっている。もし純粋に計算だけの速度を計測したい時は、mpf_out_str ()の行をコメントアウトする必要がある。
/*
The Gauss AGM algorithm using Schonhage variation.
[1] Jorg Arndt, Christoph Haenel, "Pi Unleashed", 
    pp. 93, Springer, 2001.
(C) 2003 Hironobu SUZUKI, licensed by GPL2 
*/
#include <stdio.h>
#include <gmp.h>
#include <time.h>

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

int main (int ac, char *av[])
{
  long int k, loopmax;
  mpf_t A, B, a, b, s, s_1, t, t1, t2, t3, c2;
  long int prec, dprec;
  dprec = 1000000L;		/* decimal precision */
  prec = bprec(dprec);		/* binary precision (plus alpha) */
  mpf_set_default_prec (prec);
  loopmax = 21;

  mpf_init (A);			/* big A */
  mpf_init (B);			/* big B */
  mpf_init (a);			/* a */
  mpf_init (b);			/* b */
  mpf_init (s);			/* s(n) */
  mpf_init (s_1);		/* s(n-1) */
  mpf_init (t);			/* temporary */
  mpf_init (t1);		/* temporary */
  mpf_init (t2);		/* temporary */
  mpf_init (t3);		/* temporary */
  mpf_init (c2);		/* 2 constant */
  mpf_set_ui (A, 1);
  mpf_set_ui (a, 1);
  mpf_set_ui (t1, 1);
  mpf_div_ui (B, t1, 2);
  mpf_div_ui (s, t1, 2);
  mpf_set_ui (t1, 10);
  mpf_set_ui (c2, 2);
  for (k = 1; k <= loopmax; k++)
    {
      mpf_add (t1, A, B);	/* (A+B) */
      mpf_div_ui (t, t1, 4);	/*  t = (A+B)/4 */
      mpf_sqrt (b, B);		/* b = sqrt(B) */
      mpf_add (t1, a, b);	/* (a+b) */
      mpf_div_ui (a, t1, 2);	/* a = (a+b)/2 */
      mpf_mul (A, a, a);	/* A = a * a */
      mpf_sub (t1, A, t);	/*  (A-t) */
      mpf_mul_ui (B, t1, 2);	/* B = (A - t) * 2 */
      mpf_sub (t1, B, A);	/* (B-A) */
      mpf_pow_ui (t2, c2, k);	/* 2^k */
      mpf_mul (t3, t1, t2);	/* (B-A) * 2^k  */
      mpf_set (s_1, s);
      mpf_add (s, s_1, t3);	/* s = s + (B-A) * 2^k */
    }
  mpf_add (t1, a, b);
  mpf_pow_ui (t2, t1, 2);
  mpf_mul_ui (t1, s, 2);
  mpf_div (A, t2, t1);
  mpf_out_str (stdout, 10, dprec + 10, A);
  exit(0);
}
コンパイルをする
以下に例を示す。
% cc -static -O2 -I/usr/local/include pi.c -o pi /usr/local/lib/libgmp.a
確実にpiに自分でコンパイルしたlibgmp.aをリンクさせるために、リンクするライブラリをダイレクトに指定し、ライブラリのダイナミックローディング等で違うライブラリをリンクして実行することがないように-staticをつけてpi 単独で実行できるようにしている。
実行する
以下に例を示す。
 % ./pi > pi.dat
pi.datファイルのサイズは1000014になる。e指数表現で出力されているので、ファイルは0.314159...ではじまり最後は...5628e1で終る(改行はない)。尚、先に書いたように出力中には計算誤差は含まれてはいない。尚、上記のようなプログラミングをする場合、gmpを用いて計算できる最大桁数は小数点以下約6.4億桁のはずである。それ以上の場合は、gmpで設定できる precisionのサイズを越えてしまう。また、使用メモリに関しては3300万桁(繰り返し回数26回)の計算時でメモリーを約250MB強使っていることが確認できたので、可能な最大桁数の場合、数Gバイトのレベルになると予想される。
ここまでの結論
円周率を求めるいくつかのアルゴリズムの中ではGauss AGM法が非常に高速であることがわかった。その理由は、収束が格段に良いからである。小数点以下 100万桁よりもっと出力桁数を増やす場合は、 ~~GMP標準の数値出力関数では対応できないこともわかった。~~ シェルのスタックサイズ指定を無制限にする必要がある。尚、このプログラムの出力が正しいことはBBP法を用いた別のプログラムで検証している。
further work
現在、Simon-Plouffe法を使い、円周率の計算を実行中。プログラム自体は 21.4億桁まで計算できるはずである。実行しているPentium 4 1.6GHzでは21億桁までの計算には数カ月、あるいはそれ以上かかるであろうが、その前に数千万桁のオーダーで現在のメモリの搭載量約750MBが足りなくなるはずである。
Simon-Plouffe法のプログラムはCPU/Memoryの耐久テストとして非常に良いプログラムかと思われる。

super piとの比較
ftp://pi.super-computing.org/Linux_jpからsuper_pi-jp.tar.gz (Version 2.0 of the super_pi for Linux OS) をダウンロードし、実行を行い、上記のpi.cと比較してみた。 ./pi 20とすると19回のループが実行され104万桁の計算が始まる。配布ファイルにはアルゴリズムの説明はないが、使っているアルゴリズムはGauss AGMであろう。処理速度に面白い結果が出た。Athlon 1GHz上で計算は約125秒、データの出力（10進テキスト)が約0.37秒 (合計約125.5秒)だった。一方 Athlon 上で独自コンパイルしたGMP-4.1.2をリンクするpi.cの場合、計算は約47秒とsuper-pi より、はるかに高速だった。しかし、出力に約158秒(合計約205秒)かかる。なかなか面白い結果である。この結果より円周率の計算自体はsuper-piよりもgmp-4.1.2を使った方が高速であることがわかった。
super piとの比較その2
(2006年8月3日)
gmp-4.2になってから非常に高速にpiを出力するようになった。そこでsuper piと再度比較を行なってみた。(1) にある実行コードは既にオブジェクトが古く、手もとにあるLinux2.6環境では実行できない。そこで(2) Windows 用の実行コードを入手し Pentium D 3.2GHz 上で、一方はWindows上、もう一方は VMware 上で実行することとした。結果は、gmp-4.2をリンクしたpi.cの方が ~~3.8~~ 2.4倍程高速であった。

URL
(1) ftp://pi.super-computing.org/Linux_jp
(2) http://www1.coralnet.or.jp/kusuto/PI/super_pi.html
環境
CPU Pentium D 3.2GHz
GNU/Linux 環境 VMware Workstation 5.0.1 / KNOPPIX 5.0.1 / 1GB(割り当て) / gcc-4.0.4 / gmp-4.2 / 100万桁まで計算
Windows XP 環境 Windows XP SP2 / 2GB / ライブラリ等不明 / 104万行まで計算
結果
pi.c : ~~11.64~~ 15.07 sec
super_pi : 42 sec

マルチスレッドでの計算
近年、マルチコアCPUがノートパソコンにも搭載されるようになり、マルチスレッドによる並列計算化が期待されるようになってきた。円周率計算のアルゴリズムにBBP法 (Bailey-Borwein-Plouffe 1995)を使いマルチスレッド化したプログラムを使えば、マルチコア／マルチプロセッサの能力を最大限に生かす計算ができる。しかしながらBBP計算法は計算量が多く、Gauss AGM法や Chudonovsky法のように高速には円周率を計算できない。そのため、この参考用プログラムが計算する桁数は、はるかに少ない15000桁である。円周率の並列計算例として参照して頂きたい。(2006年9月20日)
mt-bbp.c
コンパイル例
 % gcc -O3 mt-bbp.c -lgmp.a -lpthread -o mt-bbp

実行例
 % mt-bbp > pi.dat  <- デフォルトはスレッド数1
 % mt-bbp single > pi.dat   <- スレッド数1
 % mt-bbp dual   > pi.dat   <- スレッド数2
 % mt-bbp quad   > pi.dat   <- スレッド数4
 % mt-bbp max    > pi.dat   <- スレッド数16 (最大)
マシンによる比較

各種マシンによる実行時間の違いを計測してみた。

最新の情報は円周率300万桁計算時間記録表を参照してください。
% time ./pi > /dev/null
実際の数値

任意の10桁を表示する

円周率小数点以下100,000,000桁までの範囲で、任意の桁からの数字10個を表示します。

小数点以下の桁数を入力:

小数点以下の位取りについての説明

円周率は3.1415926535...となりますが、小数点以下第一位を指定したい時は、ボックスへの入力は"1"、出力は1415926535となります。小数点以下第二位を指定したい時は2となり、出力は4159265358となります。

平成16年5月3日作成
1億桁にアップ平成18年12月26日

PIベンチマーク

300万桁まで計算させてみよう

GNU/Linux上で自動的にライブラリと本体プログラムをコンパイルし、円周率 300万桁を計算させ時間を計るベンチマークのαバージョンを作りました。クライアント側である円周率を計算するプログラムは、私の作ったpi.cではなく GMPサイトにあるgmp-chudnovsky.cを使っています。またその記録をサーバへ登録する機能も出来ています。サーバ側はクライアントからのデータを受けて、簡単な集計表を作成するまでは出来ています。実行環境としてKnoppixを利用するのがお勧めですが、おおよそ標準的なGNU/Linux 環境であれば動くはずです。でもダメな時は諦めてください。(2006年8 月14日初出/2006年8月29日最終改)
gmp-4.2.2 にアップデートしました。gcc-4.2.2との組み合わせで実験したところ、さらに実行速度があがっているようです。(2007年12月2日追加)
gmp-4.3.0 にアップデートしました。AMD64 / Intel64 とも高速化しています。(2009年4月21日追加)
ダウンロード: quick start pi benchmark

ベンチマークに使う円周率計算のプログラムをgmpサイトに用意されている gmp-chudnovsky.c を使うように変更しました。Chudonovsky法(1987)を使いよく考えられたアルゴリズムでプログラムされているので非常に高速です。Intel Core Duo T2300 1.66GHzのノートパソコンですら100万桁の計算は数秒です。それもあって、これからは円周率300万桁を計算した時間を計測するようになります。それに伴いpi quick startはバージョン2に上がっています。(2006年8月18日)

Gregory and David Chudnovsky http://www.pbs.org/wgbh/nova/sciencenow/3210/04.html

結果をみてみよう

いままでサーバに送られた実行時間計測データを集計した円周率300万桁計算時間記録表を見てみましょう(2006年9月20日以前のものはクリアされています)。

マルチプロセッサ/マルチコアでのPIベンチマーク
[20071129] 最近ではマルチコアが当たり前になってきました。以前、bbpアルゴリズムを使ってマルチスレッドによる円周率を作りましたが、これを使ってプロセッサ数の違いによって円周率の計算速度がどれだけ違うかを試してみるテストをPIベンチマークに加えてみました。bbpアルゴリズムは計算量が多く、あまり高速に計算できません。しかしながら、マルチプロセッサによる計算速度の違いは計測できます。もしかすると将来、効率のよい並列に計算できるアルゴリズムを実装した時に、どれくらいの速度の違いが出るか予想する基礎データとして使えるかもしれませんし、あるいは、CPUのマルチスレッド能力を計測する指針として使えるかもしれません。当面のデータが欲しいので、協力者を募ります。実行環境はDebian GNU/Linux, ubuntuで確認していますが、他のディストリビューションで動くかどうかはわかりません。ある程度データが収集できたら結果をリストしたいと思います。
Step 1: ダウンロード

Step 2: 通常のPIベンチマークを行う。

Step 3: 通常のPIベンチマークが終了後、mt-start.shを実行、さらに mt-send.shを実行しデータをサーバに送る。
一日か二日経つといつのまにかこちらのページに反映されているはずです。

[20090421] 最新のgmp-4.3.0にアップデート、これに伴いPatchを削除
[20080702] Jason Worth Martin's Intel Core2 Patchを追加
[20071129] マルチプロセッサ/マルチコアでのPIベンチマーク追加
[20071203] マルチプロセッサ/マルチコアでのPIベンチマークランキングページ追加
なぜここのベンチマークはamdの方がintelより速いのかの考察を掲載 http://h2np.net/pi/amd_vs_intel.html (2008/02/27追加)

平成15年10月3日より
すずきひろのぶ
hironobu at h2np dot net
$Id: index.html,v 1.32 2009/04/21 04:42:50 hironobu Exp hironobu $