求平方根倒数的算法

作者:

下面这个求 \( 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.

类似文章:
编程 » C++, 算法
一个短小、高效的 C++函数,用来判断指定日期是星期几:
2007 年,我们 讨论过一个算法问题 , perfect shuffle ,据称是个微软面试题:
注: 这个问题来自 China Theory Week 2008 的 Open Problems Session。
最大回撤是一个重要的风险指标。对于对冲基金和数量化策略交易,这个指标比波动率还重要。
Google 新推出了 图片搜索 ,可直接上传图片(或者用图片链接)搜索网络上的相似图片, 例子 。估计还没多少人意识到,这玩意儿是人肉搜索的大杀器,以后大家还是少上传私人照片到公开网络。
姚期智教授给清华大学新入学研究生开的一门课,课程内容:
注: 这个问题来自 China Theory Week 2008 的 Open Problems Session。