trifle

技術メモ

Fast inverse square root

高速に  \dfrac{1}{\sqrt{x}} の近似値を計算するアルゴリズム。 Quake III Arena というビデオゲームのCGプログラムがきっかけで誕生したらしい。
めっちゃ単純で、与えられた x(32bit浮動小数点とする)に対して

  1. 浮動小数点 x を対応する 32bit のビット列 i にする。
  2. i = 0x5f3759df - ( i >> 1 );
  3. i を浮動小数点に戻す。

たったこれだけで、  \dfrac{1}{\sqrt{x}} の良い近似が得られる。

ビデオゲームにあった元のソースコードは以下のような感じらしい1

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

// 1st iteration の部分はニュートン法で精度を上げている部分なので、実質的にはその上の3行が本体ということになる。what the fuck?


仕組み

英語版 Wikipedia にちゃんと書いてあるけれど、日本語で分かりやすいものがなかったので、自分で書く。


32bit 浮動小数は、1bitの符号部 s , 8bitの指数部 e , 23bitの仮数部 m に対して、

x = (-1)s * 2e - 127 * 1.m

という形をしている。
例えば、x = 5 のとき、s = 0, e = 10000001 (=129), m = 01000000000000000000000 。
例外は、+0, -0, +INF, -INF, NaN, 非正規化数で、今は考えないことにする。また、平方根なのだから x > 0 で考える。


いま、

 x = 2^{e_x} (1 + m_x)

という数を 32bit のビット列に直すことを考える。ただし、 m_x 0 \leq m_x \lt 1 を満たし、かつ2進数表現。 例えば、 x=5 のとき  e_x = 2, m_x = 0.01
ここで、両辺2が底の対数を取ってみると、

 \log_2 {x} = e_x + \log_2 {(1 + m_x)}

さらに、 \log_2 {(1 + m_x)} m_x の線形に近似することを考える( m_x は 0から1 の範囲で小さいので)。

 \log_2 {(1 + m_x)} \approx m_x + \sigma

 \sigma は 0.0430357 がもっともあてはまりがいいらしい。これを使えば、

 \log_2 {x} \approx e_x + m_x + \sigma

となる。さて、元の  x に戻ると、 x の 32bit ビット列  I_x の指す値は

 I_x = L(e_x + B + m_x) = L(e_x + m_x + \sigma + B - \sigma) \approx L \log_2 {x} + L (B - \sigma)

と書ける。ただし、  L=2^{23}, B=127 とおいた。

同値変形して、

 \log_2 x \approx \dfrac{I_x}{L} - (B - \sigma)

となる。

求める  \dfrac{1}{\sqrt{x}} (=y) について、やはり log を考えると、

 \log_2 y = - \dfrac{1}{2} \log_2 x = - \dfrac{1}{2} \left(\dfrac{I_x}{L} - (B - \sigma)\right)

一方、 \dfrac{1}{\sqrt{x}} の 32bit ビット列を  I_y とすれば、  \log_2 y \approx \dfrac{I_y}{L} - (B - \sigma) なので、

 - \dfrac{1}{2} \left(\dfrac{I_x}{L} - (B - \sigma)\right) = \dfrac{I_y}{L} - (B - \sigma)

という関係式があることがわかった。これを同値変形して、

 I_y = \dfrac{3}{2} L(B - \sigma) - \dfrac{1}{2} I_x

マジックナンバー 0x5f3759df の正体は  \dfrac{3}{2} L(B - \sigma) のbit列(の16進数表現)2 だった。



アカデミアではない世界からこういう画期的なアルゴリズムが生まれるのは面白い。マジックナンバーというのもなんか神秘的だし。

FPUで平方根を実装したものの、よい初期値が与えられないためニュートン法を3回も繰り返すことになっており、時間がかかっていたが、これを使うと回数を減らせるかなあと期待している。