约翰·卡马克(John Carmack)的不寻常的快速反平方根(Quake III)

John Carmack在Quake III源代码中有一个特殊的function,它可以计算浮点的平方根,比常规(float)(1.0/sqrt(x))快4倍,包括一个奇怪的0x5f3759df常量。 请参阅下面的代码。 有人可以一行一行地解释这里究竟发生了什么,为什么这个工作比常规的实现要快得多?

 float Q_rsqrt( float number ) { long i; float x2, y; const float threehalfs = 1.5F; x2 = number * 0.5F; y = number; i = * ( long * ) &y; i = 0x5f3759df - ( i >> 1 ); y = * ( float * ) &i; y = y * ( threehalfs - ( x2 * y * y ) ); #ifndef Q3_VM #ifdef __linux__ assert( !isnan(y) ); #endif #endif return y; } 

仅供参考。 卡马克没有写出来。 泰耶·马蒂森(Terje Mathisen)和加里·塔罗利(Gary Tarolli)都为此付出了部分(而且非常适中)的信贷,并且还有一些其他的来源。

如何衍生神话常量是一个谜。

引用Gary Tarolli的话:

哪一个实际上是用整数做一个浮点运算 – 花费很长时间才能弄清楚这是怎么工作的,为什么这个工作起作用了,我不记得细节了。

由专家math家克里斯·洛蒙(Chris Lomont)开发的一个稍微好一些的常数,试图找出原algorithm的工作原理:

 float InvSqrt(float x) { float xhalf = 0.5f * x; int i = *(int*)&x; // get bits for floating value i = 0x5f375a86 - (i >> 1); // gives initial guess y0 x = *(float*)&i; // convert bits back to float x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy return x; } 

尽pipe如此,他最初的尝试是一个math上优秀的id的sqrt(几乎相同的常数),但实际上却比Gary最初开发的algorithm要差一些,尽pipe在math上更“纯”一些。 他无法解释为什么id是如此优秀的iirc。

当然,现在的结果是比使用FPU的sqrt(特别是在360 / PS3上)慢得多,因为float和int寄存器之间的交换会导致load-hit-store,而浮点单位可以做倒数平方根源于硬件。

它只是显示了随着底层硬件变化的性质,优化如何发展。

根据这篇写得好的文章 ,

代码的魔力,即使你不能跟着它,脱颖而出,因为我= 0x5f3759df – (i >> 1); 线。 简化的牛顿 – 拉夫森(Newton-Raphson)是一个以猜测开始并用迭代进行精简的近似。 利用32位x86处理器的本质,i,一个整数,最初被设置为想要采用整数转换的逆平方的浮点数的值。 然后我被设置为0x5f3759df,减去它本身右移一位。 右移降低了我的最低位,基本上减半了。

这是一个非常好的阅读。 这只是它的一小部分。

Greg HewgillIllidanS4给出了一个很好的math解释。 我会在这里总结一下那些不想过分细节的人。

除了一些例外,任何math函数都可以用多项式和来表示:

 y = f(x) 

可以完全转化为:

 y = a0 + a1*x + a2*(x^2) + a3*(x^3) + a4*(x^4) + ... 

其中a0,a1,a2,…是常数 。 问题是,对于许多函数,如平方根,对于精确值,这个和具有无限数量的成员,它不会以某个x ^ n结束。 但是,如果我们停止一些x ^ n,我们仍然会得到一些精确的结果。

所以,如果我们有:

 y = 1/sqrt(x) 

在这个特定的情况下,他们决定丢弃第二个以上的所有多项式成员,可能是因为计算速度:

 y = a0 + a1*x + [...discarded...] 

而现在的任务是计算a0和a1,以使y与精确值有最小的差别。 他们已经计算出最合适的值是:

 a0 = 0x5f375a86 a1 = -0.5 

所以,当你把这个方程式,你会得到:

 y = 0x5f375a86 - 0.5*x 

这与你在代码中看到的那一行是一样的:

 i = 0x5f375a86 - (i >> 1); 

编辑:其实这里y = 0x5f375a86 - 0.5*x是不一样的i = 0x5f375a86 - (i >> 1); 因为将浮点移位为整数不仅除以2还将指数除以2而引起一些其他的伪影,但是仍然归结为计算一些系数a0,a1,a2 …。

在这一点上,他们发现这个结果的精确度是不够的。 所以他们另外只做了一步牛顿迭代来提高结果的准确性:

 x = x * (1.5f - xhalf * x * x) 

他们可以在一个循环中做更多的迭代,每个迭代改进结果,直到满足要求的精度。 这正是CPU / FPU的工作原理! 但似乎只有一次迭代就足够了,这也是速度的祝福。 CPU / FPU根据需要执行尽可能多的迭代以达到存储结果的浮点数的精度,并且具有适用于所有情况的更一般的algorithm。


总而言之,他们所做的是:

使用(几乎)与CPU / FPU相同的algorithm,针对1 / sqrt(x)的特殊情况利用初始条件的改进,并且不精确地计算CPU / FPU将要提前停止的时间,从而获得计算速度。