|
多倍長演算ライブラリ(UltraPrecision) |
|
逆数法 |
| 最終更新日：2006/11/22 全面書直し |
●概要
乗算が高速化されれば、除算の高速化となるが、高速除算は、一般には、逆数法が用いられる。つまり、A / B は、A * (1 / B) とするのである。これは、FFT乗算がレガシ除算より、はるかに高速で、且つ逆数を高速に求められることが前提となる。
●逆数演算原理
ニュートン法で求めるのが普通。
y = x-1 - a とおくと、漸化式は、
xn+1 = 2 * xn - a * xn2
となる。幸い、漸化式に除算が含まれないので、FFT乗算と、加減算で実現できる。2 * xn は、レガシ乗算(多倍長*Integer なので、演算回数は、N)で行い、a * xn2 を、二回のFFT乗算で行う。
実際にニュートン法で逆数を求めて確認した。
●ニュートン漸化式
ニュートン法の原理は専門書に依られたい。逆数を求めるニュートン法には、高次収束もあり、どの方法が適しているかを検討する。
○初期値
初期値は何でも良いと言う訳ではなく、その後収束を決定付ける大切な要素である。多倍長の場合は、一般には、CPUネイティブ除算で、逆数を求め初期値にする。.NET の場合は、Double による初期値なので、15桁程度の精度となる。一般に、初期値の精度の次数倍づつ精度が上がるとされている。
○基本式
Xn+1 = 2 * Xn - A * Xn2
である。二次収束と言うべきか。ループごとに初期値精度の倍づつ、精度が向上する。つまり、15桁であれば、30、60、120・・・などとなる。
収束を見るには、上式を変形する。
d = 1 - A * Xn
と置けば、
Xn+1 = Xn * ( 1 + d)
となり、d がいかに0に近づくかとなる。演算では、dの指数部を収束判定に使用する。
○高次収束
- 3次収束
d = 1 - A * Xn
Xn+1 = Xn * (1 + d + d2)
- 4次収束
d = 1 - A * Xn
Xn+1 = Xn * (1 + d) * (1 + d2)
- 5次収束
d = 1 - A * Xn
Xn+1 = Xn * (1 + (1 + d2) * (d + d2))
それぞれの d の指数部を収束条件に使用する。
●方式の得失
演算時間が短くなるものを選ぶことになる。高次で収束回数が減少しても、ループ当りの演算時間が増加するので、一概に良いとは言えない。それぞれの演算数は、以下の表のようになる。
| 演算 | 2次収束 | 3次収束 | 4次収束 | 5次収束 |
| 乗算 | 2 | 3 | 4 | 4 |
| 加減算 | 2 | 3 | 3 | 4 |
精度が上がれば乗算時間が支配的となるので、2次収束が演算数上では有利となる。4次と5次では乗算数が同じであるので、5次が有利と予測できる。下図はシミュレータで次数と収束回数、演算時間を予測したものである。初期値は15桁。薄緑色の部分がその桁での最小時間。
ほぼ同じ傾向で、劇的な差はない。 精度が高い部分では5次が有利と出た。予想通り、4次が最も不利となっている。但し、求める精度と収束回数は非線形(量子化状態)なので、一般に関係は複雑になる。
●初期値桁数予測
シミュレーションして見た。
前半
後半
上図は、次数が 2 の場合で、列見出しの15 は、初期値が15桁の時の基準演算時間、その右の各値が設定された初期値である。それぞれの目標精度では、初期値を64〜目標精度/2 までを算出している。薄緑色の部分は、最小値で、概ね、目標精度の1/8から1/4に納まっている。いずれも、初期値15桁の演算時間より短くなっている。
下図は、次数が5の場合の後半である。
後半
初期値が15桁の場合と異なり、次数が2の方が速いようだ。初期値桁数を大きくする方式(単一初期値方式と名付ける)では、結局、次数は 2 が速そうと言える。
●多段初期値
単一初期値方式より速くなる可能性があるので、シミュレーションして見た。
上図が結果で、各次数にて算出している。1/4単位とは、目標精度の1/4づつ初期値を漸化して行く意味である。薄緑色は最小値。UltraPrecisionで筆者のマシンに置いては、次数 2 で1/4単位の多段初期値方式が最も速い(単一初期値方式よりも速い)と予測できた。100万桁では、標準方式の 4.5 倍の速さが得られることになる。
●実測値
実際に、次数が2の場合の単一初期値と多段初期値について実測した。最初、効果があまりなく、現実は甘いものではないと思ったが、プログラムにバグがあり、初期値算出時の変数の桁数が目標精度になっていたためと分かり、これを求める初期値桁数に修正したところ、下図のように全体としてシミュレーション値よりは大きくなっているがシミュレーションと同じ傾向が得られた。これにて、シミュレーション自体も結構使えると思った。多段初期値1/4単位が一番成績が良い。
図で、初期値15桁が基準となる方式、回数は、ニュートンの収束回数で、基準以外では最終段の回数。（）内の桁数は、基準の値との一致桁数で、いずれも精度に近い値となっており、演算に問題がないことが分かる。100万桁では4倍速くなっている。
●多段初期値の逆数プログラム事例
参考までに、実測に使用したプログラムを紹介する。
X の逆数を求める
PP は目標精度、NNは単位段数（4 or 8）
Private Function CReciproM(ByRef X As UltraLong) As UltraLong
Dim i, PT(), M, CP, CI As Integer
MakeMITable(PP, NN, PT, M)
'多段初期値テーブル生成
Dim Y As New UltraLong(1 / X.ToDouble) '15桁の初期値
Select Case M
Case 0
Return
CReciproE(Y, X)
'多段なし
Case Else
For i = M To 1 Step
-1 '多段処理
CP = PT(i)
TmpPrecision
= CP
Y =
CReciproE(Y, X.PartialValue(0, CP \ 8))
'n段目の初期値を求める
Y =
Y.PartialValue(0, CP \ 8)
RestoreTmpPrecision()
Next
Return
CReciproE(Y, X)
'目標精度で求める
End Select
End Function
'逆数基本関数
X0 は初期値、A の逆数を求める
Private Function CReciproE(ByRef X0 As UltraLong, ByRef A As UltraLong)
As UltraLong
Dim i As Integer
Dim d As New UltraLong(0)
Dim Xn As New UltraLong(X0)
For i = 0 To 200
d = Subt(cOne, Mul(A, Xn))
If d.IsZero OrElse d.Exponent < SeriesLimit Then
Exit For
Xn = Mul(Xn, Add(cOne, d))
Next
CVN = i
Return Xn
End Function
'多段初期値処理テーブル生成
Friend Sub MakeMITable(ByVal PP As Integer, ByVal NN As Integer, ByRef
PT() As Integer, ByRef M As Integer)
Dim C, P As Integer
C = 1
M = 0
P = PP
ReDim PT(0)
PT(0) = PP
Do
P = P \ NN
If P >= 96 Then
ReDim Preserve PT(C)
PT(C) = P
C = C + 1
M = M + 1
Else
Exit Do
End If
Loop
End Sub