倒数平方根快速算法

Posted 2020-12-13 oasisyang

倒数平方根速算法

平方根倒数速算法(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-BiasBias=2^{k-1}-1M=1+ff 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\\sigmaapprox 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)={ m 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