快速开平方根算法

人们很早就在Quake3源代码中发现了类似如下的C代码，它可以快速的求1/sqrt(x)，在3D图形向量计算方面应用很广
float invSqrt(float x)
{
float xhalf = 0.5 * x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f3759df - (i >> 1); // gives initial guess
x = *(float*)&i; // convert bits back to float
x = x * (1.5 - xhalf * x * x); // Newton step
return x;
}
在分析这段代码之前，先看看传统方法是怎么求一个数的平方根的倒数的，一般采用牛顿迭代法，为描述方便，假设输入数为a，显然需要满足a>0，sqrt是C语言的求平方根函数，为方便起见，下面用sqrt(x)的形式代替x^(1/2)
求1/sqrt(a)，用迭代法，即求方程f(x)=x^(-2)-a在f(x)=0时的解，选择适当的初始值x0，代入迭代式：
x=x-f(x)/ f'(x)
化简此式得：
x=3x/2-ax^3/2
这实际上就是上面函数倒数第二行，从函数注释也可以直接看出，这一步就是牛顿迭代，一般选择一个合适的初始值开始迭代后，迭代次数越多越接近解，换句话说就是精度越高，误差越小，当误差小于可接受值，即可获得近似结果了，就这个问题而言，初始值的选择一般要在区间(0, sqrt(3/a))，证明从略
而invSqrt这个函数厉害的地方就在于，在正式迭代开始前的三行计算已经得到了一个非常接近于解的数，因此只需一次迭代，即可得到近似值，经测试，对于常用的浮点数范围，invSqrt(x)与标准解1/sqrt(x)的最大相对误差为1.75‰，平均相对误差为0.95‰，这个精度在很多时候已经满足QuakeIII的基本要求1‰了，而invSqrt(x)的速度则比直接计算1/sqrt(x)快4倍，这对于Quake和CS这类游戏的性能是非常重要的，而且如果需要更高的精度，则将迭代那一行再重复一次就可以了，相对误差会降到百万分之一的级别，只不过速度会慢一些，接下来我们来分析下这三行代码的原理
最令人迷惑的是gives initial guess这行，即对i做了移位和减法运算，不过熟悉C语言的人应该能看出来，这个算法和float浮点数的内部表示有关，分析应该从这里入手
正式开始前先轻松下，讲些历史故事，人们在QuakeIII源码发现了这个函数，于是很自然的认为这是卡马克（John Carmack）的杰作，其中0x5f3759df这个数被称为卡马克密码，我们在下面称这个数为magic，Beyond3D.com的Ryszard Sommefeldt一直在想到底是哪个家伙写了这些神奇的代码，于是就开始找作者，John Carmack在邮件回复中明确表示不是他，也不是Michael。Terje Mathisen说他写过类似的高效代码，但上面的不是。后来猜测这个来自于一些早期黑客的算法笔记，作者究竟是谁自然也难以追查了，可以肯定的是这个家伙对计算机和高数知识都有较好理解，很聪明
2003年普渡大学的数学家Chris Lomont写了一篇文章对这段代码进行了分析。论文是英文的，地址在：