Quake 平方根倒数：准确度

Question

"magic" 计算平方根倒数的方法，显然可以追溯到 Quake 游戏，在许多资料中都有描述。维基百科上有一篇很好的文章：https://en.wikipedia.org/wiki/Fast_inverse_square_root

我特别发现以下是对算法的非常好的描述和分析：https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf

我试图在本文中复制其中一些结果，但准确性存在问题。用 C 编写的算法如下：

#include <math.h>
#include <stdio.h>

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));
  // y = y * (threehalfs - (x2 * y * y));
  return y;
}

paper 声明对于所有正正态浮点数，相对误差最多为 0.0017522874。（有关代码和第 1.4 节中的讨论，请参见附录 2。）

然而，当我 "plug-in" 数字 1.4569335e-2F 时，我得到的错误大于这个预测的公差：

int main ()
{

  float f = 1.4569335e-2F;

  double tolerance = 0.0017522874;
  double actual    = 1.0 / sqrt(f);
  float  magic     = Q_rsqrt(f);
  double err       = fabs (sqrt(f) * (double) magic - 1);

  printf("Input    : %a\n", f);
  printf("Actual   : %a\n", actual);
  printf("Magic    : %a\n", magic);
  printf("Err      : %a\n", err);
  printf("Tolerance: %a\n", tolerance);
  printf("Passes   : %d\n", err <= tolerance);

  return 0;
}

输出为：

Input    : 0x1.dd687p-7
Actual   : 0x1.091cc953ea828p+3
Magic    : 0x1.08a5dcp+3
Err      : 0x1.cb5b716b7b6p-10
Tolerance: 0x1.cb5a044e0581p-10
Passes   : 0

所以，这个特定的输入似乎违反了那篇论文中的声明。

我想知道这是论文本身的问题，还是我的编码有误。如果有任何反馈，我将不胜感激！

Answer 1

您使用了错误的幻数。

0x5f3759df是最初在Quake III中使用的值，但后来发现0x5f375a86给出了更好的结果。如果您查看您引用的论文第 40 页上的图 6.1，您会发现它使用的是改进后的常量。

这是我使用 0x5f375a86:

获得的结果

Input    : 0x1.dd687p-7
Actual   : 0x1.091cc953ea828p+3
Magic    : 0x1.08a5fap+3
Err      : 0x1.cae79153f2cp-10
Tolerance: 0x1.cb5a044e0581p-10
Passes   : 1

Answer 2

让我们尝试一小段代码来重新计算相对误差的界限，并显示它比第 561 行 thesis of Matthew Robertson. Indeed, as noticed first in the answer of @squeamishossifrage and noted in the thesis of Matthew Robertson, this implementation is the one that was disclosed in the source of Quake III. In particular, the original value of the constant of Quake III can be found in the source of Quake III, in file q_math.c 中的略大。

首先，需要调整代码以在 64 位平台上运行。唯一可能需要修改的是整数类型：long 不是独立于平台的。在我的 linux 计算机上，sizeof(long) returns 8...如第 49 页的论文中更新的那样，类型 uint32_t 将确保整数的类型是与 float.

大小相同

这是代码，由 gcc main.c -o main -lm -Wall 编译，运行由 ./main 编译：

#include <math.h>
#include <stdio.h>
#include <inttypes.h>

float Q_rsqrt(float number) {
    uint32_t i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y = number;
    i = *(uint32_t *) &y;
    i = 0x5f3759df - (i >> 1); //  0x5f3759df 0x5f375a86
    y = *(float *) &i;
    y = y * (threehalfs - (x2 * y * y));
    // y = y * (threehalfs - (x2 * y * y));
    return y;
}

int main ()
{

    printf("%ld %ld\n",sizeof(long),sizeof(uint32_t));

    uint32_t i;
    float y;
    double e, max = 0.0;
    float maxval=0;
    for(i = 0x0000000; i < 0x6f800000; i++) {
        y = *(float *) &i;
        if(y>1e-30){
            e = fabs(sqrt((double)y)*(double)Q_rsqrt(y) - 1);
            if(e > max){
                max = e;
                maxval=y;
            }
        }
    }
    printf("On value %2.8g == %a\n", maxval, maxval);
    printf("The bound is %2.12g == %a\n", max, max);

    return 0;
}

对于绑定，我获得了0.0017523386721 == 0x1.cb5d752717ep-10。如您所见，它比论文中报告的略大 (0.001752287)。使用 float 而不是 double 评估错误不会对结果产生太大影响。

Quake 平方根倒数：准确度

Quake inverse-square root: accuracy

c

floating-point

floating-accuracy

square-root