【学习笔记】快速平方根倒数算法

很以前久，看《DOOM启示录》的时候就看到，当年卡马克大神在《雷神之锤》中使用了一个神奇的数字，能够通过位操作快速计算平方根倒数$y=\frac{1}{\sqrt{x}}$。但是当时并没有深究，这两天偶然看到了这篇文章，终于将我这个多年的“未解之谜”解开了，特此做下笔记，图片取自原文。

代码实现

首先是这个神奇算法的代码：

float fast_inverse_sqrt(float x)
{
    float half_x = 0.5 * x;
    int i = *((int *)&x); // 以整数方式读取X
    i = 0x5f3759df - (i>>1); // 神奇的步骤
    x = *((float *)&i); // 再以浮点方式读取i
    x = x*(1.5 - (half_x * x * x)); // 牛顿迭代一遍提高精度
    return x;
} 

其中最最重要的就是第二步，通过这个神奇的步骤，能够得到x平方根倒数的近似值。前后穿插着用整型读取浮点数等等，实在是不知所云。但实际上，其中包含着很有意思的数学背景（每次一提到数学，我就觉得好方）。

背后的数学

首先来重新温习一下机组里面的浮点数基础知识（机组60分默默飘过）。

单精度浮点数

单精度浮点数有32bit，包括1位符号位s，8位指数位e，23位尾数位m（赞一下CSDN的markdown编辑器的latex数学公式的功能~)：

$\qquad s \space \underbrace{e\cdots e}_8\space \underbrace{m\cdots m}_{23}$

其表示的数值为：$(1+m')2^{e'}$，这里的m’代表尾数部分的数值，m’只包含小数位，即$m'=0.\underbrace{m\cdots m}_{23 }$，e’则是指数部分的数值，只有整数位，$e'=\underbrace{e\cdots e}_8-B$。为了表示负的指数，指数部分要减去一个偏移量B。
再规定M表示以整数方式解释尾数二进制位的值，E表示相同的方式解释指数部分的值。可以得到以下的转换关系：

$\qquad m'=\frac{M}{L}\qquad L=2^{23}$
$\qquad e'=E-B\qquad B=127$

那么对于一个浮点数的二进制位$s \space \underbrace{e\cdots e}_8\space \underbrace{m\cdots m}_{23}$以整数方式解读的话，其值就为：$I=EL+M$（因为平方根输入只能为正，所以默认s为0）。

平方根倒数近似计算

下面开始进入正题：
对于输入x，我们要计算的是$y=\frac{1}{\sqrt{x}}$，将该式子两边取对数：$\log_2(y)=-\frac{1}{2}\log_2(x)$。然后我们将x，y都替换为浮点数表示形式：

$\qquad \log_2((1+m'_y)2^{e'_y})=-\frac{1}{2}\log_2((1+m'_x)2^{e'_x})$
$\qquad \log_2(1+m'_y)+e'_y=-\frac{1}{2}(\log_2(1+m'_x)+{e'_x})$

观察发现，上式等号两边均有$\log_2(1+v)$形式的式子，v的范围是0到1，而该式在这个定义域内的函数图像非常接近一条直线：

因此我们近似的用$x+\sigma$替换得到：

$\qquad m'_y+\sigma+e'_y=-\frac{1}{2}(m'_x+\sigma+{e'_x})$

$\sigma$为预先计算的一个常数。接下来将数值转为整形方式读取二进制位的形式：

$\qquad \frac{M_y}{L}+\sigma+E_y-B=-\frac{1}{2}(\frac{M_x}{L}+\sigma+E_x-B)$

移项合并后得到：

$\qquad \frac{3}{2}L(\sigma-B)+M_y+LE_y=-\frac{1}{2}(M_x+LE_x)$

恰好我们又发现，在等式的两边都含有以整数方式解释浮点数的表示$I=M+LE$，进一步化简有：

$\qquad I_y=-\frac{1}{2}I_x + \kappa\qquad \kappa=\frac{3}{2}L(B-\sigma)$

由此，我们得到了一个比较有用的近似公式，$\kappa$便是神奇的数字0x5f3759df 。该公式表明计算x的平方根倒数时，我们只需要用整型方式读取x，再代入上式计算得到结果的整形解释值，然后再用浮点数方式读取即可！我当时就震惊了。

牛顿迭代提高精度

经过上一步，我们有了初步的近似结果，为了进一步提高精度，我们再用牛顿迭代法计算一次。牛顿迭代法描述的是，给定连续函数f(x)，对于方程f(x)=0，确定任意初始的$x_0$，我们可以通过下面这个迭代式计算得到其解：

$\qquad x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$

我们要计算的式子为$y=\frac{1}{\sqrt{x_0}}，x_0$为已知（这个$x_0$是指输入，不是牛顿迭代的初始值），将其转换为关于y的方程就是$f(y)=\frac{1}{y^2}-x_0=0$，根据牛顿迭代可得：

$\qquad y_{n+1}=y_n-\frac{\frac{1}{y_n^2}-x_0}{\frac{-2}{y_n^3}}$
$\qquad y_{n+1}=\frac{3}{2}y_n-\frac{1}{2}x_0y_n^3$

也就是最后一步的代码了。

更多讨论

上述推导过程中可以看出，不同次方根在我们的计算中仅仅是作为一个因子而存在的，因此将该近似算法完全可以推广到其他次方根的情况下。另外，对于64位浮点数，我们也可以采用完全相同的推导过程来得到相似的结果。这只是一个近似算法，它与真实计算值之间的偏差与神奇数字的取值有很微妙的关系，这篇博客对这个问题进行了一些讨论。

真是神奇的算法！