平方根倒數快速算法

平方根倒數速算法

平方根倒數速算法(Fast inverse square root)，經常和一個十六進制的常量 0x5f3759df聯系起來。該算法大概由上個世紀90年代的硅圖公司開發出來，后來出現在John Carmark的Quake III Arena的源碼中。

源碼：

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

准備工作

IEEE浮點數標准

IEEE浮點標准采用

\[V=(-1)^{s}×M×2^{E} \]

的形式表示一個浮點數,s是符號位,M是尾數,E是階碼.

以32位float為例子:

對於規范化值,有:

\[E=Exp-Bias\\ Bias=2^{k-1}-1\\ M=1+f\\ f \in [0,1) \]

那么對於一個浮點數x,將其各段的值按整數解釋,則有(此處默認s=0):

\[I=Exp×2^{23}+f×2^{23} \]

記:

\[L=2^{23} \\F=f×2^{23} \]

則有:

\[I=Exp×L+F \]

倒數平方根快速算法

對於函數:

\[y=\frac{1}{\sqrt x} \]

兩邊取對數,並帶入浮點數表示:

\[\log ((1+f_{y})*2^{E_y})=-\frac{1}{2}\log((1+f_{x})*2^{E_x})\\ \Longrightarrow \log(1+f_{y})+E_y=-\frac{1}{2}[\log(1+f_{x})+E_x] \]

注意到f的范圍,近似處理有:

\[\log(1+f)=\sigma +f\\ \sigma\approx 0.0430357 \]

代入化簡:

\[f_y+\sigma+E_y=-\frac{1}{2}[f_x+\sigma+E_x]\\ \Longrightarrow \frac{F_y}{L}+\sigma+Exp_y-Bias=-\frac{1}{2}[\frac{F_x}{L}+\sigma +Exp_x-Bias]\\ \Longrightarrow \frac{3}{2}L(\sigma-Bias)+F_y+L*Exp_y=-\frac{1}{2}(F_x+L*Exp_x) \]

記:

\[Bias=B\\ \zeta =\frac{3}{2}L(B-\sigma)={\rm 0x5f3759df}\\ \]

則有:

\[I_y=\zeta -\frac{1}{2}I_x \]

最后將其按浮點數編碼即可.

牛頓迭代法

利用如下的迭代式可以得到很精確的解:

\[x_{n+1}=x_{n}-\frac{f(x_n)}{f'(x_n)} \]

對於上述的計算,引入函數

\[f(y)=\frac{1}{y^2}-x_0 \]

計算有:

\[y_{n+1}=y_n(\frac{3}{2}-\frac{1}{2}x_0*y_n^2) \]

Java版本與64位版本

public static float fastFloatInverseSqrt(float x) {
        float xHalf = 0.5f * x;
        int reEncode = Float.floatToIntBits(x);
        reEncode = 0x5f3759df - (reEncode >> 1);
        x = Float.intBitsToFloat(reEncode);
        x *= (1.5f - xHalf * x * x);
        return x;
    }

public static double fastDoubleInverseSqrt(double x) {
        double xHalf = 0.5d * x;
        long reEncode = Double.doubleToLongBits(x);
        reEncode = 0x5fe6ec85e7de30daL - (reEncode >> 1);
        x = Double.longBitsToDouble(reEncode);
        x *= (1.5d - xHalf * x * x);
        return x;
    }

double fastDoubleInverseSqrt(double x){
    double xhalf=0.5 * x;
    long reEncode=*((long*)&x);
    reEncode=0x5fe6ec85e7de30da-(reEncode>>1);
    x=*((double*)&reEncode);
    x*=(1.5f-xhalf*x*x);
    return x;
}

Magic Number: 0x0x5f3759df與0x5fe6ec85e7de30da