[转自http://yueweitang.org/blog/posts/inverse-square-root-algorithm-analysis.html]

下面这个求1/\sqrt{x}的函数号称比直接调用sqrt库函数快4倍,来自游戏Quake III的源代码。

float InvSqrt(float x){
    
float xhalf = 0.5f * x;
    
int i = *(int*)&x;
    i 
= 0x5f3759df - (i >> 1);
    x 
= *(float*)&i;
    x 
= x * (1.5f - xhalf * x * x);
    
return x;
}


我们这里分析一下它的原理(指程序的正确性,而不是解释为何快)。

分析程序之前,我们必须解释一下float数据在计算机里的表示方式。一般而言,一个float数据x共32个bit,和int数据一样。其中前23位为有效数字M_x,后面接着一个8位数据E_x表示指数,最后一位表示符号,由于这里被开方的数总是大于0,所以我们暂不考虑最后一个符号位。此时

x=1.M_x 2^{E_x-127}

如果我们把计算机内的浮点数x看做一个整数I_x,那么

I_x = 2^{23}E_x+M_x

现在开始逐步分析函数。这个函数的主体有四个语句,分别的功能是:

int i = *(int*)&x; 这条语句把x转成i=I_x

i = 0×5f3759df - (i>>1); 这条语句从I_x计算I_{1/\sqrt{x}}

y = *(float*)&i; 这条语句将I_{1/\sqrt{x}}转换为1/\sqrt{x}

y = y*(1.5f - xhalf*y*y); 这时候的y是近似解;此步就是经典的牛顿迭代法。迭代次数越多越准确。

关键是第二步 i = 0x5f3759df - (i>>1); 这条语句从I_x计算I_{1/\sqrt{x}},原理:

y=1/\sqrt{x},用x=(1+m_x)2^{e_x}y=(1+m_y)2^{e_y}带入之后两边取对数,再利用近似表示\log_2(1+z)\sim z+\delta,算一算就得到

I_y = \frac{2}{3}(127-\delta)2^{23}-I_x/2

若取\delta=0.0450465679168701171875\frac{2}{3}(127-\delta)2^{23}就是程序里所用的常量0x5f3759df。至于为何选择这个\delta,则应该是曲线拟合实验的结果。

2003年Chris Lomont还写了一篇文章对这些代码进行了分析,有兴趣的读者可以下载阅读。