高速逆平方根

高速逆平方根は、任意の正の実数の平方根の逆数の近似値を高速に求めるアルゴリズム。

概要

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;
}

3DCGを描画しようとすると、計算時に扱いやすくするため、ベクトルの方向を保ったまま長さを1にする(これを正規化と呼ぶ)必要がある。この際に逆平方根を計算する必要がある。当時のCPUは平方根の計算が重かったため、精度を多少犠牲にしてでも高速化する必要があった。

3次元ベクトル $\overrightarrow{v}=(x,y,z)$を正規化するには、各成分を長さで割る(長さの逆数を掛ける)。このベクトルの長さの逆数$\frac{1}{|\mathbf{𝐯}|}$ は、三平方の定理により$\frac{1}{\sqrt{x^2+y^2+z^2}}$で計算できる。ここで逆平方根を用いる。

例

例として、1/√2=0.70710678118を計算する。

    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y = number;

ここでは、変数i、x2、y、定数1.5の宣言と定義を行っている。また、x2に2*0.5=1を、yに2を代入している。

    i = * ( long * ) &y;
    i = 0x5f3759df - ( i >> 1 );
    y = * ( float * ) &i;

yのビット列をそのままiにコピー。iには1073741824が代入される。

次に、0x5f3759dfこと1597463007-535870912=1060592025を得る。

iのビット列をyにコピー。yには0.7178746461868286が代入される。この時点で誤差は1.52%である。

    y = y * (threehalfs - ( x2 * y * y ) );

ニュートン法を適用し、さらに精度を上げる。 0.7178746461868286*(1.5-(1*0.7178746461868286*0.7178746461868286))=0.70685957213と、最終的な誤差は0.03%である。

原理

$\frac{1}{\sqrt{x}}=2^{-\frac{1}{2}{\log }_{2}x}$より、平方根の計算を対数の計算に置き換える。

この節はまだ執筆途中です。加筆、訂正して下さる協力者を募集中!