DESIGN GARDEN
NumberStorm
NumberStorm

円周率のTシャツをどうぞ.

[NEW] 最大CPU数を16まで対応できるようにしました。[20100502]

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

はじめに

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

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

準備

GMP (GNU Multi-Precision)ライブラリを使うが、ディストリビューションの デフォルトでついてくるライブラリは遅い。正確に言うと、使っているハード ウェアに最適化していない汎用の環境を前提にライブラリはコンパイルされて いるため、さらに最適化することができる。
  1. GMPのソースコードを入手する ( http://gmplib.org/ )
  2. ソースコード展開後、コンフィグレーション、コンパイル、インストールを行う
デフォルトは/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で確認していますが、他のディ ストリビューションで動くかどうかはわかりません。ある程度データが収集で きたら結果をリストしたいと思います。
    一日か二日経つといつのまにかこちらのペー ジに反映されているはずです。

  • [20100502] mt-bbp.c 最大CPU数16に対応
  • [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追加)

    CLISPで書いてみた

    BBP方式を用いた円周率計算のプログラムをClisp (ANSI Common Lisp)で書いてみた。 BBP方式は式は簡単だが計算量の多いので、速度を考えて計算精度は小数点以下100桁程度にしている。
     (setf (EXT:LONG-FLOAT-DIGITS) 330)
     (defun pi () 
      (defun s(k) 
      (* (/ 1.0L0 (expt 16.0L0 k))
       (- (/ 4.0L0 (+ 1.0L0(* 8.0L0 k)))
        (+ (/ 2.0L0 (+ 4.0L0 (* 8.0L0 k)))
         (+  (/ 1.0L0 (+ 5.0L0 (* 8.0L0 k)))
          (/ 1.0L0 (+ 6.0L0 (* 8.0L0 k))))))))
      (defun m(n)
        (if (> n 90L0)  0.0L0
            (+ (m (+ n 1.0L0))
             (s n))))
      (m 0.0L0))
     (pi)
    
    出力はこんな感じ
    [3]> 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821477L0
    
    (2011/09/07追加)
    平成15年10月3日より
    すずきひろのぶ
    hironobu at h2np dot net
    $Id: index.html,v 1.37 2011/09/07 05:13:58 hironobu Exp $