FFTとAGMによる円周率計算プログラム

[English]

詳細

これは，円周率を巨大な桁数で計算するパッケージです． これを作ったきっかけは，ある研究で作った FFT ベースの多倍長計算 ルーチンのベンチマークを行ったことです． 計算速度は Super_PI ver 1.1 @東大金田研究室 と比較して約 2.5 倍高速です． 多倍長基本演算ルーチンは四則演算と平方根です． プログラムの変更で円周率以外の計算(sqrt(2) など)もできると思います．

配布

パッケージ内のファイル (pi_fftc6)

Makefile : make file - for 32bit CPU: digits <= 6,000,000,000
Makefile_64bit : make file - for 64bit CPU
Makefile_quad : make file - for 128bit FPU
config.h : machine depending macros
fftsgx.c : FFT Package in C - split-radix - use work areas
fftsg_hx.c : FFT Package in C - split-radix - use no work areas
pi_fftca.c : PI Calc. Program - standard version - use fft*gx.c"
pi_fftcs.c : PI Calc. Program - mem save version - use "fft*g_hx.c"
pi_fftcw.c : PI Calc. Program - mem swap version - use "fft*g_hx.c"
dgt_div.c : data converter - from pi.dat to Super_PI format
README.TXT : readme file
win32bin/ : DIR
win32bin/Makefile : make file - for intel C++ compiler
win32bin/pi_cs.exe : Win32 exec. - mem save version
win32bin/pi_cs_thread.exe : Win32 exec. - use FFT level threads
win32bin/pi_cw.exe : Win32 exec. - mem swap version
win32bin/pi_cw_thread.exe : Win32 exec. - use FFT level threads
win32bin/dgt_div.exe : data converter - pi.dat => pi.txt a printout format

コンパイル方法 (pi_fftc6)

ファイル "config.h" 中にあるマクロをチェックし 必要ならば修正してください．とくに，FFT_ERROR_MARGIN は 重要なパラメータです．もし FFT_ERROR_MARGIN が非常に小さいならば 計算効率は悪くなります．もし FFT_ERROR_MARGIN >= 0.5 と設定すると 計算を間違えるかもしれません．具体的なトラブルは多倍長乗算を間違えて ニュートン反復に失敗します．FFT_ERROR_MARGIN が 10 程度でも 運がよければ計算できます．この場合， 必ず検証のための検算を行ってください．

パフォーマンス (pi_fftc5 旧バージョン)

浮動小数点演算回数 (pi_fftc5)

---- 計算桁数と FFT の長さとの関係 ----
ndigit = nfft*log_10(R), R >= 10000 or 1000
R は多倍長計算の進数です． R は DBL_ERROR_MARGIN と FFT，計算機，コンパイラの精度に依存します．

メモリー使用 (pi_fftc5)

FFT アルゴリズム

よく知られているように N 桁の乗算計算は筆算の方法でプログラムすると N^2 回の基本演算が必要になり，非常に計算が遅くなります．しかし， FFT を利用して計算すると N*log(N) 回ですみます． これは具体的には 1GHz の 1 CPU で10億桁の乗算を行った場合， 前者が数十年かかるのに対して，後者は数分で済みます． FFT による a と b の乗算手順は以下のとおりです．

  1. a，b の各桁を数列 a_j，b_j とおいて，それに対する 離散フーリエ変換 A_k，B_k を FFT で計算する．
  2. C_k = A_k * B_k を計算して C_k に対して逆 FFT を行って 畳み込み c_j を得る．
  3. c_j を c_j < 計算の進数 となるように正規化して 乗算 c を得る．
ただし，FFT による畳み込みは巡回畳み込みなので，巡回しないように ゼロをつめるなどの処理が必要です．

FFT の乗算計算の高速化に関してはさまざまなテクニックがあります． 詳しくは以下の文献等を参照してください．

AGM アルゴリズム

円周率を計算するアルゴリズムは，大きく分けて二通りの方法があります． ひとつは，無限級数を計算する方法で，有名なものにアークタンジェント公式が あります．もうひとつは，漸化式による反復公式で，AGM (算術幾何平均法) が 有名です．

ここで使っている AGM アルゴリズムは，ガウスの算術幾何平均による 完全楕円積分の計算アルゴリズムと，楕円積分に関するルジャンドルの関係式 を組み合わせる方法 (ガウスルジャンドル公式) をさらに改良したものです． 詳しい導出および証明は 2004年数理解析研究所公開講座資料(pdf 253KB) を見て下さい． 

  ---- a formula based on the AGM (Arithmetic-Geometric Mean) ----
    c = sqrt(0.125);
    a = 1 + 3 * c;
    b = sqrt(a);
    e = b - 0.625;
    b = 2 * b;
    c = e - c;
    a = a + e;
    npow = 4;
    do {
        npow = 2 * npow;
        e = (a + b) / 2;
        b = sqrt(a * b);
        e = e - b;
        b = 2 * b;
        c = c - e;
        a = e + b;
    } while (e > SQRT_SQRT_EPSILON);
    e = e * e / 4;
    a = a + b;
    pi = (a * a - e - e / 2) / (a * c - e) / npow;
  ---- modification ----
    This is a modified version of Gauss-Legendre formula
    (by T.Ooura). It is faster than original version.

この公式は，二次収束し一回の反復で桁数が倍になり，わずか 26回のループの反復で14億桁以上の結果が得られます．

参考文献

  1. E.Salamin, Computation of PI Using Arithmetic-Geometric Mean, Mathematics of Computation, Vol.30 1976.
  2. R.P.Brent, Fast Multiple-Precision Evaluation of Elementary Functions, J. ACM 23 1976.
  3. D.Takahasi, Y.Kanada, Calculation of PI to 51.5 Billion Decimal Digits on Distributed Memoriy Parallel Processors, Transactions of Information Processing Society of Japan, Vol.39 No.7 1998.
  4. T.Ooura, Improvement of the PI Calculation Algorithm and Implementation of Fast Multiple-Precision Computation, Information Processing Society of Japan SIG Notes, 98-HPC-74, 1998.

ライセンス

Links

Main Page / FFT Page