|
数学定数算出 |
|
e |
| 最終更新日：2006/06/12 修正 |
●概要
級数で求めるが、収束速度を速める工夫を行う。
●級数
eX = ∑(Xk/k!) [k = 0 to ∞]
で、X = 1 とすれば良い。そのままでは、100万桁に達するには、約21万回の演算が必要となる。
●1/N する
e = (e(1/N))N
であるから、e(1/N) を求め、N乗すれば e となるはず。今、例えば、N = 108 とすれば、e0.00000001 を計算する。これで、収束が速くなる。一般的に表せば、
d = ∑((10-m)k/k!)
とすれば、e = d^10m となる。
下表は、m の値による100万桁までの収束回数を求めたものである。級数の演算時間は期待できる。
| 0 | -4 | -8 | -16 |
| 215,151 | 121,232 | 83,955 | 51,707 |
T正直 ＞ T1/N + T冪乗 であれば、この方式は成立。
●実際の演算
d = 1 + 10-m/1 + 10-2m/(1・2) + 10-3m/(1・2・3) + ・・・・・・・・・・・・
k番目の一般項は、(10-(k-1)m・10-m)/((k-1)!・k) となる。良く見ると、赤字の部分は、一つ前の項目の値である。つまり、前の部分値の10-m を乗算し、それを k で割れば良いことになる。
ところで、(10-m)k の計算は、正直に冪乗する必要はなく、浮動小数点の指数部のみの加算で済む。そうすれば乗算はなくなり、除算も DivI であるから速い。実際には、10-m で乗算しても内部では指数部の加減算のみで行っている。しかし、オーバヘッドが多いので、直接指数部の加減算をした方が速い。
演算カーネルを示すと、
m は、10m のm
SM = 1 '累積
S = 1 '部分値
For k = 1 To 十分大きな数
S.指数部 = S.指数部 - m
S = S/k
If S.指数部 < -(目標精度) Then End
SM = SM + S
Next
SM = (SM)^10m
演算時間として考慮すべき部分は、上カーネルの赤字部分である。当然、10m乗は二進展開法で算出する。最後の冪乗による精度の劣化が気になるところであるが。
●時間予測
今、N = 10m で m を、4、8、12、16 として予測してみる。
○級数演算時間
下図は、ユーティリティのカーネル演算時間予測ツールで予測した結果である。実測値は多分これより大きくなるが目安にはなる。
変数オーダ 0 が、そのまま級数で求めた場合である。これらの時間には、冪乗算時間 T冪乗 は含まれていない。次に、T冪乗 を予測する。100万桁ではまともに行うより、4倍くらい高速化されそう。
○冪乗算時間
まともに、X10^m を行えば途方もない時間が掛かってしまう。当然、二進展開法で行う。正確に予測するために、先ず10m を二進数で表し、ビット長と 1 のビット数の合計を算出する。この数が乗算回数となる。以下の表のようになった。
| 値 | 二進表現 | 乗算回数 |
| 104 | 11000011010100000 | 23 |
| 108 | 111011100110101100101000000000 | 43 |
| 1012 | 10010001100001001110011100101010000000000000 | 58 |
| 1016 | 101100011010001010111100001011101100010100000000000000000 | 77 |
次に、カーネル演算時間予測ツールで冪乗算を予測する。
桁数の大きい部分では、級数時間に加算しても影響は小さく、T正直 ＞ T1/N + T冪乗
が成立するのでので、本法は有効と言えそう。
●実測
下図は上記の方法で実際に確認した結果である。
級数直接は、∑(1/k!) で求めたもの、10^m は、1/10^m した場合である。冪乗算の時間も含まれる。結果は、収束回数/時間/真値との一致桁数となっている。桁数が大きいところでは 10^16 の効果が大きい。さすがに、一致桁数は減少するが、予め公称精度を余計にしておけば済む。