求平方根倒数的算法

作者: , 共 1500 字 , 共阅读 0

下面这个求$ 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.098
编程 » C++
C++的浮点数转整数有四种方法,直接类型转换、round、floor、ceil。其效果如下表:
我之前一直对 Delta ($ \Delta$ )和 Gamma ($ \Gamma$ )等 Greeks 指标理解得比较模糊,今晚上用笔认真推导了一下,以下是总结。数学公式永远是最清晰的表达方式。
姚期智教授给清华大学新入学研究生开的一门课,课程内容:
注: 这个问题来自China Theory Week 2008的 Open Problems Session。