第一篇：卡马克的求平方根的倒数的程序(快速平方根倒数算法)
卡马克的求平方根的倒数的程序（快速平方根倒数算法）
2009-09-18 10:45:01| 分类： 离奇的code |字号大中小 订阅
这个程序，大概在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的浮点表示格式直接看出long类型的，然后进行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)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)(i>>1);
x=*(float *)&i;其他地方都比较简单，这三句的意思，应当是把x的浮点表示格式直接看出long类型的，然后进行0x5f3759df(i>>1);对于0x5f3759df如果看成浮点数的话，二进制位为.10111110.***11011111,其指数位为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)= 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好了。但算法思想是通用的。大家可以尝试推算一下相应的平方根算法。
点击咨询 