Quake III中的一段小代码

以下是最近热议的新开源的Quake III源码中的一段tricky code，我稍作修改后如下：

 

float CarmackInvSqrt(float x, unsigned long magicNumber) { float xhalf = 0.5f*x; int i = *(int*)&x; // get bits for floating VALUE i = magicNumber - (i>>1); // gives initial guess y0 x = *(float*)&i; // convert bits BACK to float x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy return x; }

对于其中魔法数magic number，代码作者Carmack的版本是0x5f375a86，Lomont 的版本是0x5f3759df。

发现那篇中文分析文章做得不够严谨。虽然仔细分析这个东西不是像人们想象那么神秘深奥，但也不能草草了事。在做以下之前还是要先膜拜一下作者Carmack的睿智，倒不是仅仅因为这一段代码，而是因为能写出n多的代码。

 

首先对于最后一段，尽管注释和分析都说是牛顿法，或者牛顿迭代法（Newton Iteration ），但鄙人考虑来考虑去还是感觉只能从泰勒展开（Taylor Expansion ）推导出来（两者或许有内在联系，但智识有限，目前还是没能找到迭代法的推导方法）。具体方法都是本科高数内容，如下：

 

泰勒展开式：

对y = f(x)，可以用级数逼近: y = f(a) + f'(a) (x - a) / 1! + f''(a) (x - a) / 2! + ...
其中a是在f(x)定义域上的任意实数。

在这个例子中f(x) = 1 / sqrt(x) = x^-1/2 ，于是取前二项最后结果是：

y = a^-1/2 (1.5 - 0.5·x·a^-1 )

对比代码中对应行，可以发现‘xhalf’显然对应这里的0.5·x，而‘x’则对应a^-1/2

此时再往上看，这将触及这段代码最关键部位。

这 里‘x’显然不是最初的输入变量'x'了，它经过一个处理，而这个处理正如代码和分析文章所述，是一个预估值。只是分析文章中认为是迭代法的预估值，而基 于这边讨论则应该是关于泰勒展开式的中心值a。其实根据理论，无论是迭代法预估值还是展开式的预估值，其目标都是尽可能接近输入变量'x'，以期最大限度 地接近目标，这样后续的逼近或估计都能在一个尽可能好的线性环境下进行。

 

继续上述讨论，我们实际上希望得到一个尽可能接近x^-1/2 的数，其实这也就是这个函数要返回的值。这里实际涉及的是对浮点数的处理，估计经常接触浮点优化应用场景的专业人员（个人只能想出涉及浮点的信号处理、图形游戏两个场景）已经烂熟。单精度浮点数（32bit，即这里float）的IEEE标准如下，不是很复杂：
http://en.wikipedia.org/wiki/Single_precision_floating-point_format

 

在这一步，大致已经可以了解，这个算法实际是分两步走，一是利用浮点数特性用加减和移位估计这个最终结果，然后在用展开式进行修正（同样做惯浮点优化的同学肯定对此了如指掌）。

 

至于这个估计，大致说明一下，其中有些细节（例如Lomont的暴力求精）个人感觉不是特别重要，有兴趣可以一试。

首先从浮点数原理可知，其次高8位代表对2的指数加127，这里取magic number主要要解决这个127的问题。因为这里要做的1/sqrt(x)，除法对应对数运算中的减法（就是magic number后的减号），而开根，亦即取0.5次方对应对数运算中的除法，由于恰好是除以2，所以就对应计算中的向右移位。但是为了补偿这里的127的偏 移，实际上要在被减数中体现一个127/2次方，也就是大约2⁶³ 量级，也就是这个magic number的数量级，而这个具体的数值则要根据实验再做改进（而且结果和在特定应用场景中输入变量'x'的取值范围有相当关系）。需要注意的是计算的移位和减法操作都是近似的，因为它显然会影响到后续23位小数，就运算本身来说由于本来就是近似，所以没有问题。但真是因为这种影响的所具有的随机的破坏性，导致了Carmack和Lomont两个较优解并不是两个相邻或接近的数。

 

总之，根据以上分析，这不是一个特别难理解的算法，只是如不常接触浮点优化一开始看会比较费解。当然要写好这些代码，掌握什么时候用这些方法还是需要相当深的功力的。