旭穹の陋室

快速平方根倒数

浮点数的平方根倒数1x\frac{1}{\sqrt{x}}x​1​在计算规范化向量时会被用到，规范化向量的运算在众多模拟场景中应用广泛，比如在计算机图形学中（例如，用来确定入射角和反射角以模拟光照效果）。

直接计算一个数的平方根倒数的操作非常慢，因为这需要先计算出这个数的平方根，然后再用1除以这个平方根。尽管这两步操作都是由硬件来完成的，但因为它们本身的计算过程就非常慢，所以总的计算过程也就慢。

有一个利用了浮点数在内存中存储的方式的近似算法，效果出乎意料的好。实际上，这种算法已经被实现在硬件中，因此它对于软件工程师来说并不再那么重要，但我们仍然准备详细介绍这个算法，因为它本身具有令人欣赏的优雅性和极高的教育价值。

除了方法本身，该方法的发明历史也很有趣。快速平方根倒数的计算方法被广泛的归功于游戏工作室"id Software"，该工作室在其1999年的标志性游戏"Quake III Arena"中使用了这种方法。然而，这种方法似乎是通过一种"我向一个人学习，然后那个人又向另一个人学习"的方式传播开来的，最初的源头是William Kahan（就是IEEE 754标准和Kahan求和算法的发明者）。

在2005年左右，该方法在游戏开发社区中变得流行，因为当时id Software公布了游戏的源代码。这里是源代码中的相关部分，包括注释：

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;
}

我们将一步一步地经历它所做的事情，但首先，我们需要走一小段弯路。

在电脑（或者至少是价格合理的计算器）成为日常生活的一部分之前，人们是通过查对数表来进行乘法和相关运算的——他们查找aaa和bbb的对数，将它们相加，然后找结果的反对数。

在计算1x\frac{1}{\sqrt{x}}x​1​时你也可以使用相同的技巧：

快速平方根倒数基于该恒等式，所以需要快速计算xxx的对数。事实证明，这可以通过直接将一个32位的浮点数当作整数来解释和处理。

回想一下，浮点数按顺序存储符号位（对于正值等于零），指数部分（exponent）exe_xex​和尾数部分（mantissa）mxm_xmx​，对应于：

因此其对数为：

由于mx∈[0,1)m_x \in [0,1)mx​∈[0,1)，右边的对数部分可以近似为：

近似值在区间的两端都是精确的，但为了纠正区间中段的近似偏差，我们需要通过一个小的常数值σσσ来移动或调整这个近似。因此有：

现在，记住这些近似估计，并定义L=223L=2^{23}L=223（float类型浮点数中尾数的位数）和B=127B=127B=127（指数偏移量），当我们将一个浮点数xxx的位模式重新解释为整数IxI_xIx​时，我们最终得到：

其中乘上一个整数L=223L = 2^{23}L=223等价于左移23位。

如果你调整σσσ以最小化均方误差，那么你将能得到一个异常精确的近似结果。

将浮点数x重新解释为整数（蓝色）vs 对x的对数进行缩放和平移（灰色）

现在，用近似值来表示对数，我们得到：

泰酷辣，我们现在在哪里？哦，对了，我们想计算平方根的倒数。

利用恒等式log⁡2y=−12log⁡2x\log_2y = -\frac{1}{2}\log_2xlog2​y=−21​log2​x计算y=1xy = \frac{1}{\sqrt{x}}y=x​1​，我们可以将其代入我们的近似公式，并得到：

求解IyI_yIy​有：

事实证明，我们根本不需要计算对数：上面的公式只是一个常数减去xxx的整型表示的一半。写成代码为：

i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );

我们在第一行将y重新解释为一个整数，然后在第二行将其代入上述公式，公式的第一项是魔数 $\frac{3}{2}L(B-\sigma)=$0x5F3759DF，第二项用二进制移位而不是除法进行计算。

牛顿法迭代

接下来，我们将利用公式f(y)=1y2−xf(y) = \frac{1}{y^2} - xf(y)=y21​−x和一个非常好的初始值做几次自己编码的牛顿法迭代。更新规则是：

写成代码为：

x2 = number * 0.5F;
y  = y * ( threehalfs - ( x2 * y * y ) );

初始的近似计算结果已经非常好了，所以对于游戏开发者而言，只需要进行一次迭代就足够了。仅在第一次迭代后，结果的精度就可以达到99.8%，而且可以进一步迭代以提高精度——这也是硬件中所做的：x86指令完成了其中的一些操作，并保证相对误差不超过1.5×2−121.5 \times 2^ {-12}1.5×2−12。

• Wikipedia article on fast inverse square root.