快速平方根倒数算法

原文

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

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

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数学公式的功能~)：

s e⋯e8 m⋯m23
其表示的数值为： (1+m′)2e′ ，这里的m’代表尾数部分的数值，m’只包含小数位，即 m′=0.m⋯m23 ，e’则是指数部分的数值，只有整数位， e′=e⋯e8−B 。为了表示负的指数，指数部分要减去一个偏移量B。
再规定M表示以整数方式解释尾数二进制位的值，E表示相同的方式解释指数部分的值。可以得到以下的转换关系：

m′=MLL=223
e′=E−BB=127
那么对于一个浮点数的二进制位 s e⋯e8 m⋯m23 以整数方式解读的话，其值就为： I=EL+M （因为平方根输入只能为正，所以默认s为0）。

平方根倒数近似计算
下面开始进入正题：
对于输入x，我们要计算的是 y=1x√ ，将该式子两边取对数： log2(y)=−12log2(x) 。然后我们将x，y都替换为浮点数表示形式：

log2((1+m′y)2e′y)=−12log2((1+m′x)2e′x)
log2(1+m′y)+e′y=−12(log2(1+m′x)+e′x)
观察发现，上式等号两边均有 log2(1+v) 形式的式子，v的范围是0到1，而该式在这个定义域内的函数图像非常接近一条直线：

这里写图片描述

因此我们近似的用 x+σ 替换得到：

m′y+σ+e′y=−12(m′x+σ+e′x)
σ 为预先计算的一个常数。接下来将数值转为整形方式读取二进制位的形式：

MyL+σ+Ey−B=−12(MxL+σ+Ex−B)
移项合并后得到：

32L(σ−B)+My+LEy=−12(Mx+LEx)
恰好我们又发现，在等式的两边都含有以整数方式解释浮点数的表示 I=M+LE ，进一步化简有：

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

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

xn+1=xn−f(xn)f′(xn)
我们要计算的式子为 y=1x0√，x0 为已知（这个 x0 是指输入，不是牛顿迭代的初始值），将其转换为关于y的方程就是 f(y)=1y2−x0=0 ，根据牛顿迭代可得：

yn+1=yn−1y2n−x0−2y3n
yn+1=32yn−12x0y3n
也就是最后一步的代码了。

这里写图片描述

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

真是神奇的算法！