快速平方根倒数

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

$\hat{v} = \frac{\vec{v}}{\sqrt{v_x^2+v_y^2+v_z^2}}$

直接计算一个数的平方根倒数的操作非常慢，因为这需要先计算出这个数的平方根，然后再用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;
}

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

近似对数

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

$a \times b = 10^{\log a + \log b} = \log^{-1}(\log a + \log b)$

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

$\log \frac{1}{\sqrt{x}} = -\frac{1}{2} \log x$

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

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

$x = 2^{e_x}\cdot (1 + m_x)$

因此其对数为：

$\log_2x = e_x + \log_2(1+m_x)$

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

$\log_2(1+m_x) \approx m_x$

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

$\log_2x = e_x + \log_2(1+m_x) \approx e_x + m_x + \sigma$

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

$\begin{align} I_x &= L \cdot (e_x + B + m_x) \\ &= L \cdot (e_x + m_x + \sigma + B - \sigma) \\ &\approx L \cdot \log_2(x) + L \cdot(B-\sigma) \end{align}$

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

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

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

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

$\log_2x = \frac{I_x}{L}-(B - \sigma)$

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

近似结果

利用恒等式 $\log_2y = -\frac{1}{2}\log_2x$ 计算 $y = \frac{1}{\sqrt{x}}$ ，我们可以将其代入我们的近似公式，并得到：

$\frac{I_y}{L}-(B - \sigma) \approx -\frac{1}{2}(\frac{I_x}{L}-(B - \sigma))$

求解 $I_y$ 有：

$I_y \approx \frac{3}{2}L(B-\sigma)-\frac{1}{2}I_x$

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

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

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

牛顿法迭代

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

$f'(y) = -\frac{2}{y^3} \implies y_{i+1} =y_i(\frac{3}{2} - \frac{x}{2}y_i^2)=\frac{y_i(3-xy_i^2)}{2}$

写成代码为：

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

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

拓展阅读

Wikipedia article on fast inverse square root.