Quake-III Arena (雷神之锤3)是90年代的经典游戏之一。最近,QUAKE的开发商ID SOFTWARE遵守GPL协议,公开了QUAKE-III的原代码,让世人有幸目睹Carmack传奇的3D引擎的原码。我们知道,越底层的函数,调用越频繁。3D引擎归根到底还是数学运算。那么找到最底层的数学运算函数(game/code/q_math.c),必然是精心编写的。里面有很多有趣的函数,很多都令人惊奇,估计我们几年时间都学不完。
在/code/game/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 f**k?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
函数返回1/sqrt(x),这个函数在图像处理中比sqrt(x)更有用。
注意到这个函数只用了一次叠代!(其实就是根本没用叠代,直接运算)。编译,实验,这个函数不仅工作的很好,而且比标准的sqrt()函数快4倍!要知道,编译器自带的函数,可是经过严格仔细的汇编优化的啊!
这个简洁的函数,最核心,也是最让人费解的,就是标注了“what the f**k?”的一句:i = 0x5f3759df - ( i >> 1 );
再加上y = y * ( threehalfs - ( x2 * y * y ) );
两句话就完成了开方运算!而且注意到,核心那句是定点移位运算,速度极快!特别在很多没有乘法指令的RISC结构CPU上,这样做是极其高效的。
算法的原理其实不复杂,就是牛顿迭代法,用x-f(x)/f'(x)来不断的逼近f(x)=a的根。
简单来说比如求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入x-f(x)/f'(x)后有(x+a/x)/2,现在我们选a=5,选一个猜测值比如2,那么我们可以这么算
5/2 = 2.5;
(2.5+2)/2 = 2.25;
5/2.25 = xxx;
(2.25+xxx)/2 = xxxx
...
这样反复迭代下去,结果必定收敛于sqrt(5),没错,一般的求平方根都是这么算的,但是卡马克(quake3作者)真正牛B的地方是他选择了一个神秘的常数0x5f3759df 来计算那个猜测值,就是我们加注释的那一行,那一行算出的值非常接近1/sqrt(n),这样我们只需要2次牛顿迭代就可以达到我们所需要的精度.参考资料:https://blog.csdn.net/acdreamers/article/details/12349561
对于Python,我们也可以通过类似方法计算平方根。
代码如下:
def q_sqrt(num):
t = num
t = 0x5f3759df - (t >> 1)
j=0
while not (t * t <= num-9e-6 and (t+1) * (t+1) > num+9e-6):
t = (num / t + t) / 2
#对于部分数字,所求结果将非常接近num但会有细微差别.
#对运算次数进行限制,避免死循环.
if j==t:
return t
else:
j=t
return t
import time,math
li=[12321,23432,23,34]
l=[2,3,4,5,6]
l1=[0o17,0x5f,0b0101]
t1=time.time()
for i in li:
print('%.5f'%q_sqrt(i),end='\t')
for i in l:
print('%.5f'%q_sqrt(i),end='\t')
for i in l1:
print('%.5f'%q_sqrt(i),end='\t')
t2=time.time()
print('\n')
for j in li:
print('%.5f'%math.sqrt(j),end='\t')
for i in l:
print('%.5f'%math.sqrt(i),end='\t')
for i in l1:
print('%.5f'%math.sqrt(i),end='\t')
t3=time.time()
print('\n',t2-t1,'\t',t3-t2)
输出:
111.00000 153.07514 4.79583 5.83095 1.41421 1.73205 2.00000 2.23607 2.44949 3.87298 9.74679 2.23607
111.00000 153.07514 4.79583 5.83095 1.41421 1.73205 2.00000 2.23607 2.44949 3.87298 9.74679 2.23607
0.0 0.0010008811950683594
可见,这种方法也能正确计算结果且效率较高。