高速逆平方根計算アルゴリズムの小詳解

f:id:laplus-knsn:20181205112113p:plain

この記事は数値計算 Advent Calendar 2018 5日目の記事です。

はじめに

数値解析アドベントカレンダーは大学時代を思い出させてくれるので懐かしい気持ちになります。そんな大学からの経歴についてと、本日お話することは下記の記事からの延長になるので、まずは一読をオススメします。

blog.kaneshin.co

数値解析への想い

さて、本題からズレますが、数値解析アドベントカレンダーが始まってすぐに、私の研究テーマ（実際は手伝い）だったBFGS法のアドベントカレンダーも紹介されていて、テンションがあがります。

qiita.com

参考文献にも担当教授がいたりと、この分野の数値解析分野ならまだまだ負けない気がします💪

ちなみに、BFGS はメモリ効率を高めるための Limited-memory BFGS というアルゴリズムの方がより実践では使われることが多いです。7年以上前に実装した、Limited-memory BFGSのC言語のソースコードを公開しているので、もしよければお使いください。

github.com

Logic of Fast Inverse Square Root

前回の記事で Goで高速逆平方根計算アルゴリズム - Fast Inverse Square Root - についてコードを紹介しましたが、今回は少しコードを読み解いていければと思います。

func FastInverseSqrt(x float32) float32 {
    const threehalfs = float32(1.5)

    x2 := x * float32(0.5)
    y := x
    i := *(*int32)(unsafe.Pointer(&y))
    i = 0x5f3759df - i>>1
    y = *(*float32)(unsafe.Pointer(&i))
    y = y * (threehalfs - (x2 * y * y))
    return y
}

数値解析の演算回数や誤差

const threehalfs = float32(1.5)

数値解析で研究を行うときは演算回数や情報落ち、桁落ちなどの誤差に気をつけなければなりません。

今回は const threehalfs = float32(1.5) として、前もって3/2 = 1.5 として除算をするかしないかをコンパイラ依存（定数でなくても文脈から除算を前もってコンパイル時に解釈してくれるコンパイラもあります）をしないようにしたり、わかりやすいようにマジックナンバーではなくて定数名をつけたりして保守性を高めます。研究だと他の人にソースコードを渡すことが多々あるので、ここはリーダブルなコードを心がけましょう。

また、誤差伝播も気をつけなければなりません。

ポインターのキャスティング

x2 := x * float32(0.5)
y := x
i := *(*int32)(unsafe.Pointer(&y))

32ビット同士の整数と単精度浮動少数点数のポインターキャスティングを行っています。単精度の浮動少数点数の定義は

IEEE 754 での binary32 の定義は次の通りである。

符号ビット: 1ビット指数部の幅: 8ビット仮数部の幅: 23ビット

です。この、それぞれの符号、指数部、仮数部の情報を活用したまま32ビットの整数型にしています。

IEEE 754 - Wikipedia

ある意味マジックナンバー 0x5f3759df

i = 0x5f3759df - i>>1
y = *(*float32)(unsafe.Pointer(&i))

Fast inverse square root - Wikipedia

ある意味マジックナンバーについては、First approximation of the resultに記載があります。

y = 1⁄√x を求解するために、logを取ってテイラー展開を行いそれを一次の精度で打ち切ります。これによって得られた近似式に対して、得られた数というのが 0x5f3759df になるのです。

ニュートン法の適用

y = y * (threehalfs - (x2 * y * y))

ここで、ニュートン法の反復を一度行っています。ここで、最大発生する計算誤差についても0.175%程度しか起きないため、そこそこの精度を保ったまま、一度の反復で逆平方根の近似値を得ることができます。これ以上の精度が欲しい場合は同じ計算を再度行うことによって、二回目の反復を実行します。

おわりに

数値解析のアルゴリズムは研究の精度を向上させることやより理解するためにも深掘りは重要になるので、なんとなくではなく、理論を知るのを怠らないようにしていければと思います。

at kaneshin

Free space for me.