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

高速逆平方根とは？
C言語のコード
検証
アルゴリズムの要点
[1] 逆平方根の計算を対数・指数の計算に置き換える
[2] 浮動小数点型の内部表現を利用した対数・指数の近似計算
[3] ニュートン法による収束で精度アップ
感想

高速逆平方根とは？

高速逆平方根(fast inverse square root)とは、平方根の逆数 $\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);
}

【結果】
f:id:licheng:20210206202714p:plain

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

アルゴリズムの要点

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

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

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

対数の性質より

$\log_2 \frac{1}{\sqrt x} = -\frac{1}{2} \log_2 x$ 　　　… (1)

よって

$\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ビット整数)として解釈すれば、それだけでほぼほぼ $\log_2 x$ の近似計算になります。逆にlong型のバイナリデータをfloat型として解釈すればほぼほぼ $2 ^ x$ の近似計算になります。(ただし、指数部のオフセットを考慮すること。また、long型を小数部23ビットの固定小数点数とみなすこと。)

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

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

[2.1] 対数の近似

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

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

すると
$\log_2 (x) = \log_2 ( 2 ^ e (1+m) ) = e + \log_2 ( 1+m )$

ここで $m \in [0 , 1)$ なので $\log_2 ( 1+m ) \approx m + \sigma$

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

以上より
$\log_2 (x) \approx e + m + \sigma$ 　　　　　… (4)

[2.2] σの最適値

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

$e(m) = \log_2 (1+m) - m = \dfrac{\log (1+m)}{\log 2} - m$ とすると、
$e'(m) = \dfrac{1}{\log 2} \cdot \dfrac{1}{1+m} - 1$

$e'(m) = 0$ のとき $m = \dfrac{1}{\log 2} - 1$

よって
$\sigma = \dfrac{1}{2} e \left( \dfrac{1}{\log 2} - 1 \right) = 0.043035666028$ 　　　… (5)

グラフにすると下図のようになる。
f:id:licheng:20210206190147p:plain

[2.3] 整数型での解釈

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

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

$\begin{eqnarray*} X &=& E L + M\\ &=& L (e + B + m)\\ &=& L ( ( e + m + \sigma ) + ( B - \sigma) )\\ &\approx& L \log_2 (x) + L (B - \sigma) \end{eqnarray*}$