|
技術解説(数学関数編) |
|
整数指数(XN) |
| 最終更新日：2007/04/26 更新 |
●概要
MegaLong の整数冪乗を求める。指数が整数の場合は、工夫ができる。
●高速化原理
指数関数法則
XN = X(n0+n1+n2+・・・ +nk) = Xn0・Xn1・Xn2・・・Xnk
但し、 N = n0 + n1 + ・・・ +nk とする。
を利用する。
●方法
XN で、なにもしないで、冪乗すれば、N - 1 回の乗算が必要となる。
今、N を2進展開すれば、
N = a0・20 + a1・21 + +ak・2k
となる。つまり、
XN = Xa0・1・Xa1・2・Xa2・4・・・・
となる。ai は、０ か １ なので、単純になる。例えば、N = 100 とすれば、
100 = (1100100)2 なので、22 + 25 + 26 = 4 + 32 + 64 となる。つまり、
X100 = X4・X32・X64
となる。ところで、機械的に Xk に Xk を掛けた数列を生成すると、これらは、
X、X2、X4、X8、X16、X32、X64
となるが、これらは、6回の乗算で全て求められる。必要なのは X4 ・ X32 ・ X64 なので、数列の中から選んで乗積すれば良い。つまり、
99回の乗算が、6 + 2 回の乗算
となった。約12倍高速化できたことになる。一般的には、この方法での乗算回数は、
log2(N) + log2(N) / 2 （ log2(N) / 2 は、平均として半分のビットがオンとしていると仮定）
と言える。N = 2m であれば、最も効果が得られるが、N = 2m - 1 が最悪となる(全ビットがオンなので)。
下表は冪数による本法の効果を予測したものである。速度比は、通常演算との回数比で、例えば、100000000乗(1億乗)では、270万倍の速度となる。また、この方法では、1億乗の演算を普通に行えるようになる。
●演算内容
XN で、X ＞ 0 で、10m でない場合。
N が負であれば、正にし、後で逆数を返す。
この例では、Overflow処理は省略している。
Function CIPower(X As MegaLong, N As Long) As MegaLong
Dim C As New MegaLong(1)
Dim Rf As Integer = 0
If N = 0 Then
'NOP
ElseIf N < 0 Then
Rf = 1
N = -N
EndIf
Dim BB As ULong = N
Dim CC As MegaLong = X.Copy
Do
If (BB And 1) = 1 Then
'最下位ビットチェック
C = FMul(C, CC)
'乗積
End If
BB = BB >> 1
'右に1ビットシフト
If BB = 0 Then Exit Do
CC = FMul(CC, CC)
'数列生成
Loop
If Rf = 1 Then
C = Recipro(C)
End If
Return C
End Function
●シミュレーション
予測ツールで、X の 10N 乗 の演算時間を予測した。100京乗も実用範囲であることが分かる。
●実測値
下表は、1.2345678987654321N を精度2048桁まで、IPower と繰り返し乗算で求めた時間である。この高速化の効果がよく分かる。
・下図は、整数冪数をMegaLong にした場合の例である。e = (1 + 10-m)^10m ( m → ∞) を試みた例である。1065 乗しており、60桁以上正しいネイピア数を得ているのが分かる。→ネイピア数2(2000桁の例)
下の2.71828････は、正確なネイピア数