● atan() の計算
πの値の計算に、逆正接関数 atan を使った展開式を利用するつもりなので、まず m * atan(1/n)を mpf_t型で計算します。
atan(1/n) = 1 / (1*n) - 1 / (3*n3) + 1 / (5*n5)
- 1 / (7*n7) + ...
展開式の解説
ソースプログラム
mpf_set_default_prec(256) で仮数部の精度を 256 bit に設定しました。 計算結果の小数点以下 78
桁を表示します。
i は 0, 1, 2, ... と変化します。 加算・減算の判断に使います。
n は 1, 3, 5, ... と変化します。 分母の第一項です。
n2 は ni, ni**3, ni**5, ... と変化します。 分母の第二項です。
たとえ 2**32 回ループしたとしても、i は偶数・奇数の判別にのみ使用しているため、オーバーフローを心配する必要はありません。
mi (m) の値を 1、ni (n) の値は 5 として実行しました。(1 * atan(1/5) の値)
実行結果：
i 0, n 1, n2 5, ta 0.2
0.2
i 1, n 3, n2 125, ta 0.00266667
0.197333333333333333333333333333333333333333333333333333333333333333333333333333
i 2, n 5, n2 3125, ta 6.4e-05
0.197397333333333333333333333333333333333333333333333333333333333333333333333333
i 3, n 7, n2 78125, ta 1.82857e-06
0.197395504761904761904761904761904761904761904761904761904761904761904761904762
i 4, n 9, n2 1.95313e+06, ta 5.68889e-08
0.197395561650793650793650793650793650793650793650793650793650793650793650793651
(中略)
i 52, n 105, n2 2.46519e+73, ta 3.86332e-76
0.197395559849880758370049765194790293447585103787852101517688940241033969978258
i 53, n 107, n2 6.16298e+74, ta 1.51644e-77
0.197395559849880758370049765194790293447585103787852101517688940241033969978243
i 54, n 109, n2 1.54074e+76, ta 5.95447e-79
0.197395559849880758370049765194790293447585103787852101517688940241033969978244
i 55, n 111, n2 3.85186e+77, ta 2.33887e-80
0.197395559849880758370049765194790293447585103787852101517688940241033969978244
● πの値
マチンの展開式 π/4 = 4 * atan(1/5) - atan(1/239) を使いました。
実際の計算式は π = 16 * atan(1/5) - 4 * atan(1/239) です。
実行時間を知るために clock() を使い、２つの方法を試しました。 第一の方法は前項と同様で、第二の方法は隣り合う 2 項をインラインで計算します。 終了の判別は 2 回に 1 回の割合になります。
どちらの所要時間も、ほぼ同じでした。 mpf_* 関数の実行時間比率が大きいため、ほかの部分で工夫しても大勢に影響ないようです。 このほか、隣り合う 2 項を通分して計算する方法も検討しましたが、mpf_* 演算量に大差がなく、変数も多く必要なので、試していません。
約 5000 桁のπを計算しました。 2 進数 → 10 進文字列変換は clock() で時間が測れないほど高速でした。 小数点以下 4928 桁目までの値は別に調べた値と合致していました。
6222247715 8915049530 9844489333 0963 ... ＜* 調べた値。 小数点以下 4901 桁目より。
6222247715 8915049530 9844489314 ＜* プログラム実行結果より抜粋。
malloc(), free() を使う場合、mpf_init* 関数を実行する前に malloc() する必要があります。
さもないと free() 時にコアダンプが起こります。 これは Cygwin + GMP 環境のみの制約かもしれません。