[ ホーム ] [ FFT乗算 ] [ 逆数除算 ]

上へ
ニュートン法逆数の確認
高次収束
初期値
演算方法

●概要

　基本逆数方式(初期値がDouble)と多段初期値方式について、その演算方法を紹介する。

●基本事項

• いずれの場合も、与えられた数値ではなく、その仮数部で演算する。指数は、負にすれば良い。こうすれば、収束判定が安定する。
• 値の正負は別途定める。

●演算方式

○基本逆数方式 (ReciproBasic(A))

　演算方式の確認のために用意されたもので、通常は使用しない。FFT乗算を用いている。

1. A の値を調べ、演算するまでもないかどうかを判定する。Overflow、10ⁿ など。
2. A を絶対値にし、仮数部の初めの15桁程度をDouble にする。
3. その逆数を初期値にする。
4. ニュートンで漸化する
5. 判定値が 0 か、その指数が、級数判定限界に達すれば終了。
6. 指数を、1 - 指数 にする。
   
   →判定値が 0 ：例えば、1/5 は、0.2 となり、それ以上は漸化できない。こんな場合は、判定値は完全に 0 になる。

コード例は、以下。

Dim C As New MegaLong(1)
If CPreRecipro(A, C) Then                      '演算の必要がある
   Dim i As Integer
   Dim X As MegaLong = A.Abs
   X.Exp = 1 　　　　　　　　　　　　　　　　　　　'9.99999 にする
   Dim AA As Double = X.ToDouble
   C = New MegaLong(1 / AA)
   Dim d As New MegaLong(0)
   For i = 0 To 100
      d = CSubt(cOne, CFMul(X, C))            '判定値
      If d.IsZero OrElse d.Exp <= m_SerLast Then Exit For
      C = CFMul(C, CAdd(d, cOne))     
   Next
   C.Sign = A.Sign
   C.Exp = 1 - A.Exp
End If
Return C                                             '符号、指数を調整した逆数

○多段初期値方式 (Recipro(A))

　1/4段の多段初期値にて逆数を求める。前処理は、基本と同じ。

1. 公称精度を調べ、1/4段で意味のある初期値群があるかどうか算出する
2. 1/4段化できない場合は基本逆数方式で求める。
3. 1/4段化された精度数列を元に、多段初期値演算を行う。全体としては、初期値を変数とする漸化式演算となる。

• まづ、通常のDouble初期値を求める。
• 精度数列を参照しながら、必要な精度の次の初期値を求める。
• 得た初期値を、新しい初期値にして上を繰り返す。
• 精度数列が終了しTら、目標制度で値を求める。

コード例は、以下。

'メインの内部関数
Friend Function CRecipro(ByRef A As MegaLong) As MegaLong
   Dim C As New MegaLong(1)
   If CPreRecipro(A, C) Then                    '演算の必要がある
      Dim i, CP, N, PT(0) As Integer
      MakeMITable(m_Precision, PT, N)
      If N = 0 Then 　　　　　　　　　　　　　　　　'ReciproBasic
         C = CReciproBasic(A)
      Else                                                '仮数部にて多段初期値処理
         Dim X As MegaLong = A.Abs
         X.Exp = 1                                      '9.99999 にする
         Dim AA As Double = X.ToDouble
         C = New MegaLong(1 / AA)              '開始初期値(15桁)　　　　　　　　
         For i = 0 To N - 1
            CP = PT(i)
            TmpPrecision = CP                      '一時的に中間精度に設定
            C = CReciproE(C.PartialValue(CP + 16), X.PartialValue(CP + 16))
              '数列にある精度で、変数を設定し、漸化する
            RestoreTmpPrecision()              '元の精度に戻す
         Next
         C = CReciproE(C, X)                       '目標精度で演算
         C.Sign = A.Sign
         C.Exp = 1 - A.Exp
      End If
   End If
   Return C
End Function

'1/4段精度数列生成
Private Sub MakeMITable(ByVal PP As Integer, ByRef PT() As Integer, ByRef N As Integer)
   Dim P As Integer
   N = 0
   P = PP
   Do
      P = ((P \ 4) \ 8) * 8               '8の倍数で1/4段を算出
      If P >= 64 Then
         ReDim Preserve PT(N)
         PT(N) = P
         N = N + 1
      Else
         Exit Do
      End If
   Loop
   If N > 1 Then Array.Reverse(PT)    '小さい順にする
End Sub

'ニュートン漸化式を漸化する関数
Private Function CReciproE(ByRef RA As MegaLong, ByRef A As MegaLong) As MegaLong
   Dim i As Integer
   Dim d As New MegaLong(0)
   Dim Xn As New MegaLong(RA)
   For i = 0 To 200
      d = Subt(cOne, CFMul(A, Xn))
      If d.IsZero OrElse d.Exponent <= m_SerLast Then Exit For
      Xn = CFMul(Xn, Add(cOne, d))
   Next
   Return Xn
End Function