Quake III公开源码后,有人在game/code/q_math.c里发现了这样一段代码。它的作用是将一个数开平方并取倒,经测试这段代码比(float)(1.0/sqrt(x))快4倍: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;
}
code/common/cm_trace.c中也出现了这样一段解释sqrt(x)的函数,与上面的代码唯一不同的就是这个函数返回的是number*y:/*
================
SquareRootFloat
================
*/
float SquareRootFloat(float number) {
long i;
float x, y;
const float f = 1.5F;
x = number * 0.5F;
y = number;
i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;
y = y * ( f - ( x * y * y ) );
y = y * ( f - ( x * y * y ) );
return number * y;
}
这样的代码速度肯定飞快,我就不用多说了;但算法的原理是什么呢?其实说穿了也不是很神,程序首先猜测了一个接近1/sqrt(number)的值,然后两次使用牛顿迭代法进行迭代。根号a的倒数实际上就是方程1/x^2 – a = 0的一个正实根,它的导数是-2/x^3。运用牛顿迭代公式x' = x – f(x)/f'(x),式子化简为x' = x * (1.5 – 0.5a * x^2)。迭代几次后,x的值将趋于1/sqrt(a)。
但这段代码真正牛B的是那个神秘的0x5f3759df,因为0x5f3759df – (i >> 1)出人意料地接近根号y的倒数。人们都不知道这个神秘的常数是怎么来的,只能把它当作神来膜拜。这个富有传奇色彩的常数到底咋回事,很少有人说得清楚。这篇论文比较科学地解释了这个常数。
sofa~
我来抢板凳,Matrix牛还没有回复我的问题呢,[sad],乘法口诀后的留言。
呃。0x5f375a86似乎更厉害呢。
(1.0/sqrt(x)
括号不匹配……
回复:谢谢指正,已改正
Q_rsqrt(4)=0.499154;
1.0/sqrt(4)=0.500000。不知道对精度是怎么看待的?
回ls…0.1%的精度对于quake3完全够用,毕竟只是求一个向量正交化
另外我把那篇文章翻译了一下
http://www.seraphy.cn/blog/article.asp?id=17
最后面有pdf格式下载
有兴趣去看一下吧
有趣,又迷上数学了。
另外,似乎有另一种防垃圾评论的方式,就是判断内容里的链接数量和连续的评论内容相似度。刚刚不知道是因为数学不好,还是超时而被拒绝了。[lol]
您好。看了几篇文章,很有启发!尤其是位运算的。
我有点疑问?是怎么样判断代码的执行效率的呢?几倍或是什么的?
有什么工具,或者是?
抱歉,本人新来乍到,还在为Noip奋斗着……
能否请教,针对Noip,或者根高,该怎么样提高呢?
其实卡马克也是通过暴力方法得到这个常数的
请问 这句得到的是什么 i = * ( long * ) &y;
不要告诉我这段代码是卡马克本人写的,这不是他本人写的,这不是……
对6楼提到的翻译过的论文很感兴趣,可惜链接点不开了…
可惜可惜…
昨天在研究这个。。。
牛,真是太牛了
csdn上的确有……
我测试了一下,神奇的是double下还是(float)(1.0/sqrt(i*1.0f));这个快
Function InvSqrt(x#):Dim xhalf#,i&:xhalf=x/2:i=&H5F3759DF-x2:x=i*(1.5-xhalf*x*x):InvSqrt=x:End Function
Private Sub Form_Load():MsgBox InvSqrt(4):End Sub
10duowei
真…神奇!
抱歉,本人新来乍到,还在为Noip奋斗着…
报告,文中“牛顿迭代法”的链接失效了
http://www.matrix67.com/blog/archives/361
Good!
写的太经典了,我竟无言以对。