高速逆平方根
fast inverse square root
昔のコンピュータグラフィックで有用だった
光の反射を高速に計算できる
計算コストの高い浮動小数点除算をバイパスできる
マジックナンバー0x5f3759dfが利用されていた
code:fast_invsqrt.jl
function fast_invsqrt(x::Float32)
MAGIC = 0x5f3759df
i = reinterpret(Int32, x) # 浮動小数点数を整数として解釈
i = MAGIC - (i >> 1) # マジックナンバーを使って初期値を求める
y = reinterpret(Float32, i) # 整数を浮動小数点数として解釈
y = y * (1.5f0 - 0.5f0 * x * y * y) # ニュートン法で精度を上げる
return y
end
# using Plots
# x = 1f0:0.1f0:10f0 # Float32型で定義する
マジックナンバーの導出
1. 浮動小数点数の内部表現を利用して$ \log_2{x} を得る 対数を近似する
実数$ x を指数部$ e と仮数部$ m\in[0,1) を用いて表す
$ x\coloneqq 2^e(1+m)
これにより、$ \log_2{x} は:
$ \log_2{x}=\log_2{2^e(1+m)}=e+\log_2{(1+m)}
ここで、$ m\in[0,1) であるので、$ \log_2{(1+m)} はパラメータ$ \delta を用いて表現できる:
$ \log_2{(1+m)}\approx m+\delta
ゆえに、$ \log_2{x}\approx e+m+\delta を得る
この$ \delta がマジックナンバーを生むあんも.icon
$ \delta を最適化する: $ \delta = \log_2{(1+m)} - m
$ \delta の最適値として、$ \max{|\log_2{(1+m)} - m|} の最小化問題を考える
$ m\in[0,1) であるので、$ \delta を$ f(m) = \log_2{(1+m)} - m の最大値の半分になるように定めればよい
$ f(m) = \frac{\log(1+m)}{\log2}-m
$ f'(m) = \frac{1}{\log2}\cdot\frac{1}{1+m} - 1
$ f(m) は$ m = \frac{1}{\log2} - 1 で最大値をとる
したがって、$ \delta = \frac{1}{2} f\left(\frac{1}{\log2} - 1\right) と定められる
おおよそ0.0430356660279671あんも.icon
2. int型で解釈する
Float32の構成
符号部: 1bit
指数部: 8bit
仮数部: 23bit
$ X \coloneqq \mathrm{int}(x) と定義する
指数部$ E = e+B : B=127
仮数部$ M = m\times L : L=2^{23}
$ X = EL+M = L(e+B+m) = L\{(e+m+\delta) + (B-\delta)\}
$ \approx L\log_2(x) + L(B - \delta)
ゆえに$ \log_2(x) \approx \frac{X}{L} - (B - \delta) を得る
3. マジックナンバーを導出する: $ y = \frac{1}{\sqrt{x}}
$ \log_2y = -\frac{1}{2}\log_2x
$ \iff \frac{Y}{L} - (B - \delta) = -\frac{1}{2}\left(\frac{X}{L} - (B - \delta)\right)
$ \iff Y = \frac{3}{2}L(B-\delta) - \frac{1}{2}X
$ \frac{3}{2}L(B-\delta) がマジックナンバーとなる
$ mgc = 3/2 * 2^(23) * (127 - 0.0430356660279671); UInt32(Float32(mgc))
0x5f37bc80が得られる