個人的な課題として、πの値を最速で求める方法を探しています。具体的には、M_PI
のような#define
定数を使ったり、ハードコードで数値を入力したりしない方法を考えています。
以下のプログラムは、私が知っている様々な方法をテストしています。インライン・アセンブリのバージョンは、理論的には最も速いオプションですが、明らかにポータブルではありません。他のバージョンと比較するためのベースラインとして、このバージョンを含めました。私のテストでは、ビルトインを使用した場合、GCC 4.2では 4 * atan(1)
バージョンが最速でした、これは atan(1)
を定数に自動折りたたむからです。fno-builtinを指定すると、
atan2(0, -1)` バージョンが最速になります。
以下は、メインのテストプログラム(pitimes.c
)です:
#include <math.h>
#include <stdio.h>
#include <time.h>
#define ITERS 10000000
#define TESTWITH(x) { \
diff = 0.0; \
time1 = clock(); \
for (i = 0; i < ITERS; ++i) \
diff += (x) - M_PI; \
time2 = clock(); \
printf("%s\t=> %e, time => %f\n", #x, diff, diffclock(time2, time1)); \
}
static inline double
diffclock(clock_t time1, clock_t time0)
{
return (double) (time1 - time0) / CLOCKS_PER_SEC;
}
int
main()
{
int i;
clock_t time1, time2;
double diff;
/* Warmup. The atan2 case catches GCC's atan folding (which would
* optimise the ``4 * atan(1) - M_PI'' to a no-op), if -fno-builtin
* is not used. */
TESTWITH(4 * atan(1))
TESTWITH(4 * atan2(1, 1))
#if defined(__GNUC__) && (defined(__i386__) || defined(__amd64__))
extern double fldpi();
TESTWITH(fldpi())
#endif
/* Actual tests start here. */
TESTWITH(atan2(0, -1))
TESTWITH(acos(-1))
TESTWITH(2 * asin(1))
TESTWITH(4 * atan2(1, 1))
TESTWITH(4 * atan(1))
return 0;
}
そして、x86とx64のシステムでのみ動作するインラインアセンブリのもの(fldpi.c
)です:
double
fldpi()
{
double pi;
asm("fldpi" : "=t" (pi));
return pi;
}
そして、私がテストしているすべての設定をビルドするビルドスクリプト(build.sh
)です:
#!/bin/sh
gcc -O3 -Wall -c -m32 -o fldpi-32.o fldpi.c
gcc -O3 -Wall -c -m64 -o fldpi-64.o fldpi.c
gcc -O3 -Wall -ffast-math -m32 -o pitimes1-32 pitimes.c fldpi-32.o
gcc -O3 -Wall -m32 -o pitimes2-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -fno-builtin -m32 -o pitimes3-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -ffast-math -m64 -o pitimes1-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -m64 -o pitimes2-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -fno-builtin -m64 -o pitimes3-64 pitimes.c fldpi-64.o -lm
様々なコンパイラのフラグでテストする以外に(最適化が違うので32ビットと64ビットも比較しました)、テストの順番を入れ替えてみたりもしました。しかし、やはりatan2(0, -1)
版が毎回上位を占めています。
モンテカルロ法」1は、前述のように、いくつかの素晴らしいコンセプトを適用していますが、明らかに、最速ではありませんし、長いショットでもありませんし、合理的な尺度でもありません。また、どのような精度を求めるかにもよります。私が知っている最速のπは、桁がハードコードされているものです。円周率][2]や[円周率[PDF]][3]を見ると、たくさんの公式がある。
ここでは、1回の反復で約14桁という高速に収束する方法を紹介します。現在最速のアプリケーションである[PiFast][4]は、この式をFFTと組み合わせて使っています。コードは簡単なので、式だけ書いておきます。この式は、Ramanujanがほぼ発見し、Chudnovskyが発見したものです。実際に彼が数十億桁の数を計算した方法である--だから、無視できる方法ではない。この式はすぐにオーバーフローしてしまうので、階乗を割るのであれば、そのような計算を遅らせて項を削除するのが得策である。
。
。
のところです、
。
以下、[Brent-Salaminアルゴリズム][10]です。Wikipediaによると、aとbが "close enough" のとき、(a + b)² / 4t がπの近似となるそうです。
let pi_2 iters =
let rec loop_ a b t p i =
if i = 0 then a,b,t,p
else
let a_n = (a +. b) /. 2.0
and b_n = sqrt (a*.b)
and p_n = 2.0 *. p in
let t_n = t -. (p *. (a -. a_n) *. (a -. a_n)) in
loop_ a_n b_n t_n p_n (i - 1)
in
let a,b,t,p = loop_ (1.0) (1.0 /. (sqrt 2.0)) (1.0/.4.0) (1.0) iters in
(a +. b) *. (a +. b) /. (4.0 *. t)
最後に、πゴルフ(800桁)なんていかがでしょう?160文字です!
int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}
[2]: http://functions.wolfram.com/Constants/Pi/ "円周率"; [3]: http://functions.wolfram.com/PDF/Pi.pdf "円周率計算式"; [4]: http://numbers.computation.free.fr/Constants/PiProgram/pifast.html "PiFast";
[10]: http://mathworld.wolfram.com/Brent-SalaminFormula.html "Brent Salamin Formula";
ここでは、高校時代に習った円周率を計算するテクニックの一般的な解説をします。
モンテカルロ法とは、ランダムな過程ではすぐに導き出せない答えを導き出すための統計的手法です。
正方形を描き、その中に四分円(半円の4分の1)を刻む(半径が正方形の辺と等しい四分円は、正方形をできるだけ埋めるようにする)
正方形に向かってダーツを投げ、その着地点を記録してください。もちろん、ダーツは正方形の内側に着地したが、半円の内側に着地したか?この事実を記録する。
この作業を何度も繰り返すと、半円の中にある点の数と投げられた総数の比があることがわかります。この比をxと呼びます。
正方形の面積はr倍なので、半円の面積はx倍r倍r倍(つまりx倍rの2乗)であることが推論されます。したがって、x×4でπが求まる。
これはすぐに使える方法ではない。しかし、モンテカルロ法の良い例です。また、このような方法を使えば、自分の計算能力では解決できないような多くの問題が解決できることに気づくかもしれません。