这个程序，大概在2006年看到的，当时进行了分析，主要分析了位运算那部分：如何利用浮点数的位表示法快速计算一个近似值。而对于整个迭代计算原理并没有仔细看，今天再重新分析一下，并把当时那部分分析写的更明了些。
原贴如下：----------------------------------------------------------------------------------------------------------
发信人: interma ( 4PZP | 抓紧时间 | OFG P ), 信区: C_Cpp
标 题: [zz] 源自Quake3的快速求InvSqrt()函数
发信站: 兵马俑BBS (Mon Dec 4 13:11:36 2006), 本站(202.117.1.8)
http://games.solidot.org/article.pl?sid=06/12/04/0210205&from=rss
"人们很早就在Quake3源代码中发现了如下的C代码，它可以快速的求1/sqrt(x)，在3D图形向量计算方面应用很广。
float InvSqrt(float x){
float xhalf=0.5f*x;
long i=*(long*)&x;
i=0x5f3759df - (i>>1);
x=*(float *)&i;
x=x*(1.5f-xhalf*x*x);
return x;
}
Beyond3D.com的Ryszard Sommefeldt一直在想到底是哪个家伙写了这些神奇的代码？2003年Chris Lomont还写了一篇文章（PDF)对这些代码进行了分析。毫无疑问写出这些代码的人绝对是天才。"
是John Carmack？Michael Abrash？John Carmack在邮件回复中明确表示不是他，也不是Michael，可能是Terje Matheson。于是侦探Ryszard又向Terje Mathisen寻求答案。Terje说他写过类似的高效代码，但上面的不是。他猜测可能在MIT的旧文档中可以找到。Ryszard的求证之路显然无止境。
☆──────────────────────────────────────☆
phylips 于 Tue Dec 5 10:04:42 2006 提到：
分析了一下，关于计算过程还是可以理解的
但是关于其中的数学依据，只能看出一部分
主要思想利用了浮点数的表示法，ieee的标准
其中关键部分是
long i=*(long*)&x;
i=0x5f3759df - (i>>1);
x=*(float *)&i;
其他地方都比较简单，这三句的意思，应当是把x的浮点表示格式直接看出long类型的，
然后进行0x5f3759df - (i>>1)运算，之后再直接把long类型的看成float类型的。其中
的数学道理，我大致看了下
对于浮点数的表示，在ieee里float类型是按照符号.指数.尾数格式进行的
其中符号位1位，指数位8位，尾数23位，先只考虑了指数部分，
long i=*(long*)&x是把浮点数看成了与他具有相同位格式的long型数
i=0x5f3759df - (i>>1)
对于0x5f3759df如果看成浮点数的话，其指数位为(0101)(1111)0->1011111=190
设原来的xIEEE表示的指数位为a则新的结果的指数位变成了190-a/2；注意如果要转换成
m*2^n格式，需要将指数a-127；所以得到的新结果用数学表示的指数位为190-a/2-127=
63-a/2
而x的原来如果从ieee格式变成数学表示应为a-127
比较63-a/2与a-127
可以看出指数之间的关系为63-a/2==（a-127)*(-1/2)
所以经过上述三句，x应该已经是1/sqrt(x)的近似值
然后下面的应该是继续进行了精化，具体依据还没搞懂
【 在 interma ( 4PZP | 抓紧时间 | OFG P ) 的大作中提到: 】
☆──────────────────────────────────────☆
xiaocui 于 Tue Dec 5 10:06:34 2006 提到：
厉害。
【 在 phylips (爱立佛) 的大作中提到: 】
☆──────────────────────────────────────☆
NemoK 于 Tue Dec 5 11:58:01 2006 提到：
re!
【 在 phylips (爱立佛) 的大作中提到: 】
☆──────────────────────────────────────☆
zqfalcon 于 Tue Dec 5 12:01:02 2006 提到：
呵呵大牛啊，佩服
【 在 phylips 的大作中提到: 】
---------------------------------------------------------------------------------------------------------------------
首先看牛顿迭代法的原理
若函数f（x）在点的某一临域内具有直到（n+1）阶导数，则在该邻域内f（x）的n阶泰勒公式为：
f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2! +… +...fn(x0)(x- x0)n/n!....
解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点附近展开成泰勒级数 f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2! +… 取其线性部分，作为非线性方程f(x) = 0的近似方程，即泰勒展开的前两项，则有f(x0)+f'(x0)(x－x0)=f(x)=0 设f'(x0)≠0则f(x0)+f'(x0)(x－x0)=0 其解为x1=x0－f(x0)/f'(x0) 这样，得到牛顿法的一个迭代序列：x(n+1)=x(n)－f(x(n))/f'(x(n))。
也就是说x(n+1)=x(n)－f(x(n))/f'(x(n))，是用来计算f(x)=0的根的迭代公式，如果x(n)与真根足够靠近的话，那么只需要少数几次迭代，就可以得到满意的解。
我们现在要求的是一个数的平方根的倒数，设这个数为a，那么我们要计算的就是x= a^(-0.5)。首先我们要把这个计算需求转换为f(x)=0的形式。
x = a^(-0.5)
x^(-2) = a
x^(-2)-a = 0
这样就转化成了f(x) = 0的形式，其中x就是a的平方根的倒数。
f(x) = x^(-2)-a 则f'(x) = -2*x^(-3)
根据牛顿迭代公式可以得到
x(n+1) = x(n) - (x(n)^(-2) - a)/(-2*x(n)^(-3))
x(n+1) = x(n) + (x(n)^(-2) - a)/(2*x(n)^(-3))
x(n+1) = x(n) + (x(n)/2) - (a*x(n)^3)/2
x(n+1) = 1.5*x(n) - (0.5*a*x(n)^3)
这样我们就得到了迭代计算的公式，也就是源程序中的 x=x*(1.5f-xhalf*x*x);
下面关键的一步就是计算x(n)的一个真根足够靠近的近似值。主要思想利用了浮点数的表示法，ieee754标准
其中关键部分是
long i=*(long*)&x;
i=0x5f3759df - (i>>1);
x=*(float *)&i;
其他地方都比较简单，这三句的意思，应当是把x的浮点表示格式直接看出long类型的，然后进行0x5f3759df - (i>>1)运算，之后再直接把long类型的看成float类型的。其中的数学道理，我大致看了下
对于浮点数的表示，在ieee里float类型是按照符号.指数.尾数格式进行的，其中符号位1位，指数位8位，尾数23位，先只考虑指数部分。
long i=*(long*)&x;是把浮点数看成了与他具有相同位格式的long型数，
i=0x5f3759df - (i>>1);对于0x5f3759df如果看成浮点数的话，二进制位为.10111110.01101110101100111011111,其指数位为10111110=190 。
设原来的x的IEEE表示的指数位为a则新的结果的指数位变成了190-a/2；注意如果要转换成m*2^n格式，需要将指数减去127；所以得到的新结果用数学表示的指数位为190-a/2-127= 63-a/2 ，而x的原来的指数a，如果从ieee格式变成数学表示应为a-127 。
比较63-a/2与a-127，可以看出指数之间的关系为63-a/2==（a-127)*(-1/2)
所以经过上述三句，x应该已经是1/sqrt(x)的近似值。当然上面只是从指数位进行了简单分析，如果要细化到尾数，还需要看具体效果。
------------------------------------------------------------------------------------------------------------------------------------------
以下摘自：http://tbl.javaeye.com/blog/231425 雷神之锤III》里求平方根倒数的函数（快速平方根(倒数)算法）
至于那个0x5f3759df，呃，我只能说，的确是一个超级的Magic Number。
那个Magic Number是可以推导出来的，但我并不打算在这里讨论，因为实在太繁琐了。简单来说，其原理如下：因为IEEE的浮点数中，尾数M省略了最前面的1，所以实际的尾数是1+M。如果你在大学上数学课没有打瞌睡的话，那么当你看到(1+M)^(-1/2)这样的形式时，应该会马上联想的到它的泰勒级数展开，而该展开式的第一项就是常数。下面给出简单的推导过程：
对于实数R>0，假设其在IEEE的浮点表示中， 指数为E，尾数为M，则：
R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)
将(1+M)^(-1/2)按泰勒级数展开，取第一项，得：
原式
= (1-M/2) * 2^(-E/2)
= 2^(-E/2) - (M/2) * 2^(-E/2)
如果不考虑指数的符号的话， (M/2)*2^(E/2)正是(R>>1)， 而在IEEE表示中，指数的符号只需简单地加上一个偏移即可， 而式子的前半部分刚好是个常数，所以原式可以转化为：
原式 = C - (M/2)*2^(E/2) = C - (R>>1)，其中C为常数
所以只需要解方程：
R^(-1/2)
= (1+M)^(-1/2) * 2^(-E/2)
= C - (R>>1)
求出令到相对误差最小的C值就可以了
上面的推导过程只是我个人的理解，并未得到证实。而Chris Lomont则在他的论文中详细讨论了最后那个方程的解法，并尝试在实际的机器上寻找最佳的常数C。
而Chris Lomont则在他的论文中详细讨论了最后那个方程的解法，并尝试在实际的机器上寻找最佳的常数C。
所以，所谓的Magic Number，并不是从N元宇宙的某个星系由于时空扭曲而掉到地球上的，而是几百年前就有的数学理论。只要熟悉NR和泰勒级数，你我同样有能力作出类似的优化。
在GameDev.net上有人做过测试，该函数的相对误差约为0.177585%，速度比C标准库的sqrt提高超过20%。如果增加一次迭代过程，相对误差可以降低到e-004的级数，但速度也会降到和sqrt差不多。据说在DOOM3中，Carmack通过查找表进一步优化了该算法，精度近乎完美，而且速度也比原版提高了一截。
值得注意的是，在Chris Lomont的演算中，理论上最优秀的常数（精度最高）是0x5f37642f，并且在实际测试中，如果只使用一次迭代的话，其效果也是最好的。但奇怪的是，经过两次NR后，在该常数下解的精度将降低得非常厉害（天知道是怎么回事！）。经过实际的测试，Chris Lomont认为，最优秀的常数是0x5f375a86。如果换成64位的double版本的话，算法还是一样的，而最优常数则为0x5fe6ec85e7de30da（又一个令人冒汗的Magic Number - -b）。
这个算法依赖于浮点数的内部表示和字节顺序，所以是不具移植性的。如果放到Mac上跑就会挂掉。如果想具备可移植性，还是乖乖用sqrt好了。但算法思想是通用的。大家可以尝试推算一下相应的平方根算法。
评论