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

Rushyo.. 71

仅供参考.卡马克没有写下来.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.

仅供参考.卡马克没有写下来.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.

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

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

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的所有方式将转到但更早停止,因此获得计算速度.

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

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

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

我很好奇看到常量是浮点数所以我只是编写了这段代码并用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