逛知乎的时候看到这么个问题：
其中有个答案提到了一种计算浮点数的平方根倒数的快速算法，其实我以前就曾经思考过大数平方根的计算方法，然后…就没有然后了。
神奇的地方在于这个快速算法中有个诡异的常数，据说至今没有人知道它最早是怎么来的。
代码示例
下面这段代码来源于上面那个答案。
1 | float InvSqrt(float x) |
这个算法说到头只有这么几行而已，而恰恰神奇的就是这么几行代码可以非常快速地实现平方根倒数。
整个算法的关键就是第5行的那个诡异常数以及第7行的**牛顿迭代**了。
算法原理
要理解这个算法，首先要明白牛顿迭代，而牛顿迭代法的基础是泰勒级数。
泰勒级数
定义：如果$f(x)$在点$x=x_0$处具有任意阶导数，则幂级数：
$$
\begin{align}
f(x)&=\sum^{\infty}_{n=0}\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n\
&=f(x_0)+f’(x_0)(x-x_0)+\frac{f’’(x_0)}{2!}(x-x_0)^2+…+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+…
\end{align}$$
称为$f(x)$在$x=x_0$处的泰勒级数。
泰勒展开在数学上还有很多的应用…然而，高数学过去那么多年了，有啥用我早忘了…0.0
这里只要知道泰勒级数可以用来近似计算函数值就行了。
后面部分明显是越往后越小的，收敛性应该是要有条件来严格证明的…这个我也忘了。
牛顿迭代
牛顿迭代解线性方程，就是把线性方程$f(x)=0$线性化的一种近似方法。把$f(x)$在$x_0$的某个邻域内展开成泰勒级数，取其线性部分（前两项）：
$$
f(x)=f(x_0)+f’(x_0)(x-x_0)
$$
令其为零，就可以作为原本$f(x)=0$方程的近似方程了。只要$f’(x_0)\neq0$，方程的解就可以写成：
$$
x_1=x_0-\frac{f(x_0)}{f’(x_0)}
$$
这样，就得到牛顿迭代法的迭代关系式了：
$$
x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)}
$$
理论上来说，是不是我任意取一个满足$f’(x_0)\neq0$的数，都可以把它作为第1个近似解，然后通过不断地进行牛顿迭代，直到精度达到要求即解出了比较靠谱的方程解？
回到目标算法上面来，假设已知的数是$a$，要求它的平方根的倒数，就是解方程：
$$
(\frac{1}{x})^2=a
$$
即解方程：
$$
f(x)=\frac{1}{x^2}-a=0
$$
那么：
$$
f’(x)=-\frac{2}{x^3}
$$
代回到牛顿迭代式中去：
$$
\begin{align}
x_{n+1}&=x_n-\frac{f(x_n)}{f’(x_n)}\
&=x_n-\frac{\frac{1}{x^2_n}-a}{-\frac{2}{x^3_n}}\
&=x_n+\frac{x_n-ax^3_n}{2}\
&=\frac{3}{2}x_n-\frac{a}{2}x^3_n
\end{align}
$$
这个就是上面程序第7行的那个牛顿迭代式了。上面那个程序段只迭代了一次，其实精度是不够的，一般迭代2次，跟标准值的误差就不大了。
当然，迭代次数越多，精度越高，相应地耗时也是越长。
诡异的常数
最后，还有个问题没解决……0x5f3759df这个常数是什么鬼啊？
上面已经讨论过了，通过牛顿迭代法，我们可以在第n个近似解的基础上推出第n+1个精度更高的近似解。这里唯一的漏洞就是：我们必须得到第1个近似解！这样才能够往下递推！
而上面第5行的这段代码：
1 | i = 0x5f3759df - (i >> 1); |
就是用来计算第1个近似解的。
好吧，还是不明白这个数是怎么来的…
据史料记载：
浮点数的平方根倒数常用于计算正规化矢量。3D图形程序需要使用正规化矢量来实现光照和投影效果，因此每秒都需做上百万次平方根倒数运算，而在处理坐标转换与光源的专用硬件设备出现前，这些计算都由软件完成，计算速度亦相当之慢；在1990年代这段代码开发出来之时，多数浮点数操作的速度更是远远滞后于整数操作，因而针对正规化矢量算法的优化就显得尤为重要。
1999年的《雷神之锤III竞技场》就借助了这个神奇的算法。
《雷神之锤III》的代码直到QuakeCon 2005才正式放出，但早在2002年（或2003年）时平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了。最初人们猜测是卡马克写下了这段代码，但他在询问邮件的回复中否定了这个观点，并猜测可能是先前曾帮id Software优化雷神之锤的资深汇编程序员Terje Mathisen写下了这段代码；而在Mathisen的邮件里他表示在1990年代初他只曾作过类似的实现，确切来说这段代码亦非他所作。现在所知的最早实现是由Gary Tarilli在SGI Indigo中实现的，但他亦坦承他仅对常数R的取值做了一定的改进，实际上他也不是作者。Rys Sommefeldt则在向以发明MATLAB而闻名的Cleve Moler查证后认为原始的算法是Ardent Computer公司的Greg Walsh所发明，但他也没有任何决定性的证据能证明这一点。
目前不仅该算法的原作者不明，人们也仍无法明确当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验：他编写了一个函数，以在一个范围内遍历选取R值的方式将逼近误差降到最小，以此方法他计算出了线性近似的最优R值0x5f37642f（与代码中使用的0x5f3759df相当接近），但以之代入算法计算并进行一次牛顿迭代后，所得近似值与代入0x5f3759df的结果相比精度却仍略微更低；而后Lomont将目标改为遍历选取在进行1－2次牛顿迭代后能得到最大精度的R值，并由此算出最优R值为0x5f375a86，以此值代入算法并进行牛顿迭代后所得的结果都比代入原始值（0x5f3759df）更精确，于是他的论文最后以“原始常数是以数学推导还是以反复试错的方式求得”的问题作结。在论文中Lomont亦指出64位的IEEE754浮点数（即双精度类型）所对应的魔术数字是0x5fe6ec85e7de30da，但后来的研究表明代入0x5fe6eb50c7aa19f9的结果精确度更高（McEniry得出的结果则是0x5FE6EB50C7B537AA，精度介于两者之间）。在Charles McEniry的论文中，他使用了一种类似Lomont但更复杂的方法来优化R值：他最开始使用穷举搜索法，所得结果与Lomont相同；而后他尝试用带权二分法寻找最优值，所得结果恰是代码中所使用的魔术数字0x5f3759df，因此McEniry确信这一常数或许最初便是以“在可容忍误差范围内使用二分法”的方式求得。
先复习一下计算机中**浮点数的表示法**
然后回到上面讨论过的第1个近似解，第1个近似解越接近正确答案，则显然牛顿迭代法的效果会越好。
考虑已知的浮点数$a$在计算机中的实际表示为：
$$(1+F)*2^E$$
那么它的平方根倒数为：
$$(1+F)^{-\frac{1}{2}}*2^{-\frac{E}{2}}$$
猜测一下：第5行代码中，将i右移1位，再用常数减掉，可以看成是将指数除2之后再取反。至于尾数…早都变得不知道哪去了。
所以，还是没明白这个常数是怎么来的，我只能猜测到这里了。
这里把Chris Lomont写的那篇论文也一并贴出来：FAST INVERSE SQUARE ROOT