Fast Inverse Square Rootの一般化を考えてみた

Fast Inverse Square Rootという、Quakeという昔のゲームで使われた超絶すごいアルゴリズムがありまして、計算コストが割りと高い$y = 1/\sqrt x$を超高速で計算できます。

上のビデオで詳しく解説されていますが、この手法は浮動小数点数の仕様をハックしきっていて最高なので、その一般化版を考えてみました。つまり、

$y = x^{-1/2}, x^{-1/4}, x^{-1/8}...$ に対する解です。

これ全部$1/2^n$の形なので$log_2$したときにbit shiftの形に落とし込めるのがミソです。

ビデオで紹介されている箇所はなるべく端折りますが、以下でやっていきます。

（符号少数点数の符号部は全て正とします。）

$\Gamma = y^{-1/2^n} $
$\Leftrightarrow log_2{\Gamma} = {-1/2^n} log_2{y}$
$\Leftrightarrow log_2{(1+\frac{M_{\Gamma}}{2^{23}})2^{E_{\Gamma}-127}} = {-1/2^n} log_2{(1+\frac{M_{y}}{2^{23}})2^{E_{y}-127}}$
$\Leftrightarrow 1/2^{23} (M_{\Gamma} + 2^{23} E_{\Gamma}) + \mu - 127 =
{-1/2^n} (1/2^{23} (M_{y} + 2^{23} E_{y}) + \mu - 127)$
$\Leftrightarrow M_{\Gamma} + 2^{23} E_{\Gamma} = (2^{23} + 2^{23-n})(127-\mu) -1/2^n (M_{y} + 2^{23} E_{y})$

右辺第一項は定数ですが、適当にn=5,そして動画でも出ていた$\mu=0.043$で計算した場合、
(2^23+2^18)(127-0.043) = 1098273521.66 ≒ 1098273522 = 0x417652F2
となります。
最後の数式右辺第二項の$-1/2^n$が、-してn回bit shiftに当たります。

次にNewton法のパートですが、

$y = \frac{1}{x^{1/2^n}}$

を変形して

$y^{-2^n} - x \equiv f(y)$

とします。

導関数は、

$f'(y) = -2^n y^{-2^n-1}$

です。

その為、値更新のための漸化式は、

$y_{n+1} = y_n - \frac{f(y)}{f'(y)}$
$= y_n - \frac{y_n^{-2^n} - x}{-2^n y_n^{-2^n-1}}$
$= y_n + (1 + 1/2^n) y_n - 1/2^n x y_n^{2^n + 1}$

です。

なので、全体のコードとしては以下です（n=5の場合）。

float calcPower(float x, int power){
    tmp = x // puts x^(2^n) of x;
    ans = 1;
    while (power > 0){
        if (power % 2 == 1){
            ans *= tmp;
        }
        tmp *= tmp;
        power >>= 1;
    }
    return ans;
}

float Q_rsqrt( float number )
{
 const int n = 5;
 long i;
 float x2, y, y_pow;
 const float newton1 = 1.03125; // (1 + 1/2^n)
 const float newton2 = 0.03125; // (1/2^n)
 const float newton3 = 33;      // (2^n)

 x2 = number * 0.5F;
 y  = number;
 i  = * ( long * ) &y;                       // evil floating point bit level hacking
 i  = 0x417652F2 - ( i >> n );               // what the fuck?
 y  = * ( float * ) &i;
 y_pow = calcPower(y, newton3)
 y  = y * ( newton1 - newton2 * x * y_pow  );   // 1st iteration
 y_pow = calcPower(y, newton3) // 2nd iteration, this can be removed
 y  = y * ( newton1 - newton2 * x * y_pow  );   // 2nd iteration, this can be removed

 return y;
}

工夫として、累乗の部分は高速化を加えてみました。a^5を計算するとき、aaaaaだと乗算5回ですが、a^2を計算、a^2a^2を計算、a^4aを計算だと乗算3回で済むことを利用しています。

0x417652F2のマジックナンバーや、newwton1, 2, 3の変数は各自でnに合わせて計算していただければと思います。また、これ申し訳ないのですが手元にCの環境がないため検証できません。どなたか検証していただければ嬉しいです。