求平方根倒数的算法

作者: , 共 1217 字

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

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

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

分析程序之前,我们必须解释一下 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$$

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

  1. int i = *(int*)&x; 这条语句把\( x\) 转成\( i=I_x\)
  2. i = 0x5f3759df - (i>>1); 这条语句从\( I_x\) 计算\( I_{1/\sqrt{x}}\)
  3. y = *(float*)&i; 这条语句将\( I_{1/\sqrt{x}}\) 转换为\( 1/\sqrt{x}\)
  4. 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{3}{2}(127-\delta)2^{23}-I_x/2$$

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

Q. E. D.

类似文章:
相似度: 0.115
编程 » C++
C++的浮点数转整数有四种方法,直接类型转换、round、floor、ceil。其效果如下表:
相似度: 0.079
编程 » C++, 算法
一个短小、高效的 C++函数,用来判断指定日期是星期几:
我之前一直对 Delta (\( \Delta\) )和 Gamma (\( \Gamma\) )等 Greeks 指标理解得比较模糊,今晚上用笔认真推导了一下,以下是总结。数学公式永远是最清晰的表达方式。
姚期智教授给清华大学新入学研究生开的一门课,课程内容:
注: 这个问题来自China Theory Week 2008的 Open Problems Session。