卡马克快速平方根
其实很早以前就知道卡神创造了一个神奇的方法可以在线性时间内完成求一个数平方根倒数的运算，但是没有太在意。最近因为一个事把这茬想起来了，也就顺便一口气完成这个的证明吧。
卡马克源码
快速平方根在quake3源码的game/code/q_math.c文件中，这个文件提供了随机数，碰撞检测，求法向量，角度弧度互换等众多功能，但其中最霸气的，还是这个快速求出平方根倒数的Q_rsqrt函数。(Quake-III Arena的下载地址)。
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
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}
简化代码
float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i >> 1); // 计算第一个近似根
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x); // 牛顿迭代法
return x;
}
数学推导
牛顿迭代算法
该算法的本质其实就是牛顿迭代法（Newton-Raphson Method，简称NR），而NR的基础则是泰勒级数（Taylor Series）。NR是一种求方程的近似根的方法。首先要估计一个与方程的根比较靠近的数值，然后根据公式推算下一个更加近似的数值，不断重复直到可以获得满意的精度。其公式如下：
则方程\(f(x)\)的第n+1个近似根为
NR的公式可以看出：理论上，只要迭代次数足够，总能逐渐趋近于准确解。NR最关键的地方就是估算第一个近似根。如果第一个根和结果很接近，那么只需要几次迭代就能得到精度足够的解。
求平方根倒数
求平方根的倒数，其实就是求\(\frac{1}{x^2}-a=0\)的解。将该方程用牛顿迭代法解开：
如果把1/2放到括号内，就得到倒数第二行代码x = x*(1.5f - xhalf*x*x)
接下来，就要求第一个近似根\(x_n\)，这也是整个函数最神奇的地方，因为只用一次迭代过程就得到了这个解：
i = 0x5f3759df - (i >> 1);
float数据类型
float数据类型位数作用：
| bits | 31 30 … 0 |
|---|---|
| 31 | 符号位 |
| 30-23 | 共8位，保存指数(E) |
| 22-0 | 共23位，保存尾数(M) |
所以，32位浮点数用十进制实数表示就是：
开根号取倒数就是：
0x5f3759df是哪来的
语句(i >> 1)的作用就是将指数除以2，得到\(2^{-\frac{E}{2}}\)的部分。
对于实数R>0，假设在IEEE的浮点数表示中，指数为E，尾数为M，则其原理可以简述为为IEEE的浮点数中，尾数M省略了最前面的1，所以实际的尾数是1+M。如果你在大学上数学没忘干净，那么当你看到\((1+M)^{-\frac{1}{2}}\)这样的形式时，应该会马上联想的到它的泰勒级数展开，而该展开式的第一项就是常数。
将\((1+M)^{-\frac{1}{2}}\)按泰勒级数展开，取第一项
在不考虑指数符号的情况下，\(\frac{M}{2}2^{\frac{E}
{2}}\)正是(R>>1)。
在IEEE浮点数中，指数的符号可以考位移实现(指数位的阶码 = 阶码真值 + 127)，而\(2^{-\frac{E} {2}}\)是常数项。所以原式可以化为：
之后，只要求出另误差最小的C值即可。
[转]后续
普渡大学的数学家Chris Lomont看了以后觉得有趣，决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。Lomont也是个牛人，在精心研究之后从理论上也推导出一个最佳猜测值，和卡马克的数字非常接近, 0x5f37642f。卡马克真牛，他是外星人吗？
传奇并没有在这里结束。Lomont计算出结果以后非常满意，于是拿自己计算出的起始值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确的求得平方根。结果是卡马克赢了… 谁也不知道卡马克是怎么找到这个数字的。
最后Lomont怒了，采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数字要好上那么一丁点的数字，虽然实际上这两个数字所产生的结果非常近似，这个暴力得出的数字是0x5f375a86。
Lomont为此写下一篇论文，”Fast Inverse Square Root”。
下载地址: http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf
http://www.matrix67.com/data/InvSqrt.pdf
c/c++标准库的修改
从C++09开始，STL的平方根算法，就已经更换为了更快的方法。其实自己实验就可以知道，卡马克方法迭代一次比C++09标准库块，但是迭代两次速度就是差不多的了。
