#define
個人的な挑戦として、π の値を取得する最も速い方法を探しています。具体的には、などの定数を使用したり、 の数値をハードコーディングしたりしない方法を使用していますM_PI
。
以下のプログラムは、私が知っているさまざまな方法をテストします。インライン アセンブリ バージョンは、理論的には最速のオプションですが、明らかに移植性はありません。他のバージョンと比較するための基準としてこれを含めました。私のテストでは、組み込みの場合、バージョンは4 * atan(1)
GCC 4.2 で最速です。これは、自動的に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;
}
fldpi.c
そして、 x86 および x64 システムでのみ機能するインライン アセンブリ ( ) は次のとおりです。
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
のモンテカルロ法は、前述のように、いくつかの優れた概念を適用していますが、明らかに最速ではありません。また、どのような合理的な基準で見ても、はるかに高速ではありません。また、どのような精度を求めているかによっても異なります。私が知っている最も高速なπは、数字がハードコードされているものです。円周率そして円周率[PDF]、数式がたくさんあります。
ここでは、反復ごとに約 14 桁ですぐに収束する方法を示します。ピファスト現在最も高速なアプリケーションは、FFTでこの式を使用しています。コードは簡単なので、式だけ書きます。この式は、ラマヌジャンとチュドノフスキーによって発見された実際に彼は数十億桁の数字を計算したので、無視できる方法ではありません。数式はすぐにオーバーフローし、階乗を割り算しているので、そのような計算を遅らせて項を削除するのが有利です。
どこ、
以下はブレント・サラミンアルゴリズムWikipedia には、aとbが「十分に近い」場合、 (a + b)² / 4t がπ の近似値になると書かれています。「十分に近い」とはどういう意味かわかりませんが、私のテストでは、1 回の反復で 2 桁、2 回で 7 桁、3 回で 15 桁になりました。もちろん、これは倍数の場合なので、表現によっては誤差が生じる可能性があり、実際の計算はより正確である可能性があります。
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);}