John Carmack不寻常的快速反向平方根(Quake III)

Question

John Carmack在Quake III源代码中有一个特殊的功能,它计算浮点的平方根,比常规的快4倍(float)(1.0/sqrt(x)),包括一个奇怪的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;
}

Answer 1

仅供参考.卡马克没有写下来.Terje Mathisen和Gary Tarolli都为此获得了部分(并且非常适度)的信用,并且还记入其他一些来源.

这个神话常数是如何产生的,这是一个神秘的东西.

引用加里塔罗利:

实际上这是以整数形式进行浮点计算 - 花了很长时间才弄清楚这是如何以及为什么这样做的,我不记得细节了.

由专家数学家(Chris Lomont)开发的一个稍好的常数,试图弄清楚原始算法是如何工作的:

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;
}

尽管如此,他最初尝试的数学'上级'版本的id's sqrt(几乎相同的常数)证明不如最初由加里开发的版本,尽管在数学上更"纯粹".他无法解释为什么id是如此优秀的iirc.

Answer 2

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

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

Answer 3

Greg Hewgill和IllidanS4给出了一个很好的数学解释链接.对于那些不想过多介绍细节的人,我会在这里总结一下.

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

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);因为将float移动为整数不仅除以2而且将指数除以2并导致其他一些伪影,但它仍然归结为计算一些系数a0,a1,a2 .......

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

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

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

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

使用(差不多)与CPU/FPU相同的算法,利用1/sqrt(x)的特殊情况下的初始条件的改进,并且不计算精确CPU/FPU的所有方式将转到但更早停止,因此获得计算速度.

Answer 4

根据这篇不错的文章写了一会儿...

即使你不能遵循它,代码的神奇之处在于i = 0x5f3759df - (i >> 1); 线.简化后,Newton-Raphson是一种近似值,它以猜测开始并通过迭代进行精炼.利用32位x86处理器的特性,i,一个整数,最初设置为您想要使用整数转换的反转平方的浮点数的值.然后我将其设置为0x5f3759df,减去自身向右移位一位.右移降低了i的最低位,基本上将其减半.

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

Answer 5

我很好奇看到常量是浮点数所以我只是编写了这段代码并用google搜索弹出的整数.

    long i = 0x5F3759DF;
    float* fp = (float*)&i;
    printf("(2^127)^(1/2) = %f\n", *fp);
    //Output
    //(2^127)^(1/2) = 13211836172961054720.000000

看起来常量是"2 ^ 127的平方根的整数近似,其浮点表示的十六进制形式更为人所知,0x5f3759df" https://mrob.com/pub/math/numbers-18.html

它在同一个网站上解释了整个事情.https://mrob.com/pub/math/numbers-16.html#le009_16