离奇的code

卡马克的求平方根的倒数的程序(快速平方根倒数算法)

2009年9月18日 阅读(1,029)

这个程序,大概在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好了。但算法思想是通用的。大家可以尝试推算一下相应的平方根算法。

You Might Also Like