高速逆平方根(fast inverse square root)のアルゴリズム解説

高速逆平方根とは？

高速逆平方根(fast inverse square root)とは、平方根の逆数 1x√$\frac{1}{\sqrt{x}}$ を高速に計算するアルゴリズムです。平方根の逆数は逆平方根とも呼ばれます。逆平方根はベクトルの正規化などに用いられるので、これを高速に計算できるアルゴリズムには大きなご利益があります。

参照： Fast inverse square root - Wikipedia

C言語のコード

高速逆平方根の関数を示します。0x5F3759DF という謎のマジックナンバーが目を引きます。(後述)

float invsqrt(float x)
{
    // y = 1/sqrt(x) = 2 ^ (-1/2 * log2(x))
    long X, Y;
    float y;
    X = *(long*)&x;
    Y = 0x5F3759DF - (X >> 1); // Magic number!
    y = *(float*)&Y;

    // Newton's method
    const float threehalfs = 1.5F;
    float x2 = x * 0.5F;
    y = y * (threehalfs - (x2 * y * y));   // 1st iteration
//  y = y * (threehalfs - (x2 * y * y));   // 2nd iteration

    return y;
}

検証

1から100までの数について、上記の関数で逆平方根を計算し、標準ライブラリのsqrt関数を用いた計算と比較しました。

【テストコード】

#include <stdio.h>
#include <math.h>

int main(void)
{
    double ave_error = 0;
    double max_error = 0;

    for (double x = 1; x <= 100; x += 1)
    {
        // use sqrt()
        double y = 1 / sqrt(x);

        // fast inverse square root
        float fx = (float)x;
        float fy = invsqrt(fx);

        printf("%f, %f, %f\n", x, y, fy);

        // error
        double error = fabs((fy - y) / y) * 100;
        ave_error += error;
        if (error > max_error) max_error = error;
    }
    ave_error /= 100;

    printf("ave_error = %f, max_error = %f\n", ave_error, max_error);
}

【結果】

グラフがほとんど重なっています。誤差は最大0.175%、平均0.088%ほどでした。このように、わずかな計算量にもかかわらずかなり精度の良いアルゴリズムです。

アルゴリズムの要点

高速逆平方根のアルゴリズムの要点は以下の3点です。

1. 逆平方根の計算を対数・指数の計算に置き換える
2. 浮動小数点型の内部表現を利用した対数・指数の近似計算
3. ニュートン法による収束で精度アップ

[1] 逆平方根の計算を対数・指数の計算に置き換える

対数の性質より

log21x√=−12log2x${\log }_2\frac{1}{\sqrt{x}}=-\frac{1}{2}{\log }_{2}x$　　　… (1)

よって

1x√=2−12log2x$\frac{1}{\sqrt{x}}=2^{-\frac{1}{2}{\log }_{2}x}$　　　　　　… (2)

これにより、逆平方根の計算を対数・指数の計算に置き換えることができます。

[2] 浮動小数点型の内部表現を利用した対数・指数の近似計算

C言語等で用いられる float型 (32ビット浮動小数点数) の内部表現のバイナリデータ形式に着目すると、対数・指数の簡易な近似計算ができます。

参照： 単精度浮動小数点数 - Wikipedia

float型のバイナリデータは最上位ビットから順に以下にようになっています。

• 符号ビット: 1ビット
• 指数部: 8ビット (ただし、+127のバイアスがある。)
• 仮数部: 23ビット (ただし、[0,1) ではなく [1, 2) の範囲を表す。)

指数部が上位にあることに着目してください。ここで大きなトリックがあります。float型のバイナリデータをlong型(32ビット整数)として解釈すれば、それだけでほぼほぼ log2x${\log }_{2}x$ の近似計算になります。逆にlong型のバイナリデータをfloat型として解釈すればほぼほぼ 2x$2^x$ の近似計算になります。(ただし、指数部のオフセットを考慮すること。また、long型を小数部23ビットの固定小数点数とみなすこと。)

参照： 2^x と log2(x) の高速な近似計算 - 滴了庵日録

これをもう少し詳しく見ていきます。

[2.1] 対数の近似

実数 x$x$ を、指数部 e$e$ と仮数部 m∈[0,1)$m\in [0,1)$ を用いて次のように表す。

x=2e(1+m)$x=2^e(1+m)$　　　　　… (3)

すると
log2(x)=log2(2e(1+m))=e+log2(1+m)${\log }_2(x)={\log }_2(2^e(1+m))=e+{\log }_2(1+m)$

ここで m∈[0,1)$m\in [0,1)$ なので log2(1+m)≈m+σ${\log }_2(1+m)\approx m+\sigma$

ただし σ$\sigma$ は近似を調整するパラメータであり、 σ≈0.0430357$\sigma \approx 0.0430357$ で最適となる。(後述)

以上より
log2(x)≈e+m+σ${\log }_2(x)\approx e+m+\sigma$　　　　　… (4)

[2.2] σの最適値

σ$\sigma$の最適値として、ここでは誤差の絶対値の最大値が最小となる値を考える。
m∈[0,1)$m\in [0,1)$ の範囲での log2(1+m)−m${\log }_2(1+m)-m$ の最大値の半分の値を σ$\sigma$ とすればよい。

e(m)=log2(1+m)−m=log(1+m)log2−m$e(m)={\log }_2(1+m)-m={\displaystyle \frac{\log (1+m)}{\log 2}}-m$ とすると、
e′(m)=1log2⋅11+m−1$e^{\prime }(m)={\displaystyle \frac{1}{\log 2}}\cdot {\displaystyle \frac{1}{1+m}}-1$

e′(m)=0$e^{\prime }(m)=0$ のとき m=1log2−1$m={\displaystyle \frac{1}{\log 2}}-1$

よって
σ=12e(1log2−1)=0.043035666028$\sigma ={\displaystyle \frac{1}{2}}e({\displaystyle \frac{1}{\log 2}}-1)=0.043035666028$　　　… (5)

グラフにすると下図のようになる。

[2.3] 整数型での解釈

さて、実数 x$x$ のfloat型バイナリデータをlong型として解釈した整数を X$X$とする。

指数部 E=e+B$E=e+B$ ただし B=127$B=127$ はバイアス、
仮数部 M=m×L$M=m\times L$ ただし L=223$L=2^{23}$ とすると、(3) (4) より

X===≈EL+ML(e+B+m)L((e+m+σ)+(B−σ))Llog2(x)+L(B−σ)$\begin{array}{rcl}X& =& EL+M\\ & =& L(e+B+m)\\ & =& L((e+m+\sigma )+(B-\sigma ))\\ & \approx & L{\log }_2(x)+L(B-\sigma )\end{array}$

よって
log2(x)≈XL−(B−σ)${\log }_2(x)\approx {\displaystyle \frac{X}{L}}-(B-\sigma )$　　　　　… (6)

[2.4] 逆平方根の計算とマジックナンバー0x5F3759DF

y=​1x√$y=\frac{1}{\sqrt{x}}$ とすると、(1) より

log2(y)=−12log2(x)${\log }_2(y)=-\frac{1}{2}{\log }_2(x)$

ここで、y$y$ のfloat型バイナリデータをlong型として解釈した整数を Y$Y$とすると、(6) より

YL−(B−σ)=−12(XL−(B−σ))${\displaystyle \frac{Y}{L}}-(B-\sigma )=-{\displaystyle \frac{1}{2}}({\displaystyle \frac{X}{L}}-(B-\sigma ))$

よって、
Y=32L(B−σ)−12X$Y=\frac{3}{2}L(B-\sigma )-\frac{1}{2}X$　　　　　… (7)

ここで、32L(B−σ)$\frac{3}{2}L(B-\sigma )$ = 0x5F3759DF が謎のマジックナンバーの正体である。

ただし、(5) で求めた σ$\sigma$ の値から計算すると 0x5F37BCB6 となり、一般的に知られているマジックナンバー 0x5F3759DF と微妙に食い違います。この理由はよく分かりません。

(7) で求めた Y$Y$ をfloat型のバイナリデータとして解釈すれば、y=​1x√$y=\frac{1}{\sqrt{x}}$ の近似値が求まる。
これはまさに (2) の近似計算である。

[3] ニュートン法による収束で精度アップ

[2]の方法で求めた近似値から、さらに誤差を0に収束させるためにニュートン法を用いる。

誤差が0であれば、y=​1x√$y=\frac{1}{\sqrt{x}}$ より y−2−x=0$y^{-2}-x=0$

そこで関数 f(y)=y−2−x$f(y)=y^{-2}-x$ を定義し、 f(y)=0$f(y)=0$ の解をニュートン法で求める。

まず導関数は f′(y)=−2y−3$f^{\prime }(y)=-2y^{-3}$

よって、y$y$ の近似値 yn$y_n$に対し、次の漸化式でより正確な近似値 yn+1$y_{n+1}$ が得られる。

yn+1==yn−f(yn)f′(yn)yn(32−12xy2n)$\begin{array}{rcl}{y}_{n+1}& =& y_n-{\displaystyle \frac{f(y_n)}{f^{\prime }(y_n)}}\\ & =& y_n(\frac{3}{2}-\frac{1}{2}x{y}_n^2)\end{array}$

[2]の方法で求めた近似値を初期値 y0$y_0$ として、この漸化式の計算を繰り返すことでより正確な近似値が得られます。先に示したC言語のコードでは、この漸化式を1回のみ計算しています。(2回目の計算はコメントアウトしています。)

感想

これ思いついたやつ、頭おかしい。
わずかな行数のコードに込められた数学がすごすぎる。