Quake反平方根算法：准确性

Question

“魔法”计算反平方根的方法据说可以追溯到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;
}

这篇论文指出，对于所有正常浮点数，相对误差最大为0.0017522874（请参见附录2中的代码以及第1.4节中的讨论）。

然而，当我“插入”数字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

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

我想知道这是否是论文本身存在的问题，还是我的编码有误。我会非常感激任何反馈！

- alias

---

查看论文，附录A.2，最大误差是通过抽样计算得出的。因此，如果发现误差略大于最大误差，我不会太担心... - francis

@francis 不完全是。该程序测试从0x00800000到0x7f7fffff的每个32位值，几乎涵盖了所有正浮点数的范围。 - r3mainer

我认为 @francis 是正确的；采样是问题所在。所述相对误差适用于 Quake 中使用的原始魔法数字，这也是我使用的数字。(即，问题实际上关乎原始魔法数字及其相关的相对误差；而不是针对它的改进。) - alias

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

- r3mainer

---

确实，你是对的。这个 Quake III 的原始值可以在 Quake III 的源代码中，在文件 q_math.c 的第 561 行找到。那一行的注释并不打算包含在问答网站上... 这个常量的历史和 64 位 IEEE754 双精度浮点数的最优值 (0x5fe6eb50c7b537a9) 可以在 [维基百科] 上找到，引用了 Matthew Robertson 的工作，他还报告了双精度和四重精度的值！ - francis

我不确定我是否同意。如果您查看论文的第4.7节，改进的魔数0x5f375a86的最大相对误差为0.0017512378；这与原始魔数的误差是不同的。因此，我认为我的问题仍然存在，关于原始的Q_rsqrt函数的相对误差。 - alias

@LeventErkok 我明白你的意思。从情况来看，作者使用 Q_rsqrt() 来指代原始的 Quake III 代码，而使用 rsqrt() 来指代改进版。这两个函数都没有产生第4.7节中引用的结果。也许你可以尝试联系作者以获得澄清。 - r3mainer

为了测试Q_rsqrt()函数在我的64位计算机上，我觉得最好使用unsigned int代替long。 我发现Q_rsqrt()的最大误差为0x1.cb5d752717ep-10=0.001752338672。 它是由浮点数0x1.dd678p-49获得的。 正如@LeventErkok所指出的那样，它比引用的界限更高... - francis

我试图找到作者，但是找不到他的电子邮件地址。也许有一天他会看到这个！@francis 感谢进一步的分析。如果你能把那个评论变成一个答案，我很乐意接受！ - alias

Answer 2

让我们尝试一小段代码，重新计算相对误差的界限，并显示它略大于Matthew Robertson的论文中的界限。事实上，正如@squeamishossifrage的答案中首先注意到并在Matthew Robertson的论文中指出的那样，这个实现是Quake III源代码中公开的实现。特别地，Quake III常量的原始值可以在Quake III的源代码文件q_math.c的第561行找到。

首先，需要调整代码以在64位平台上运行。唯一可能需要修改的是整数类型：long不是平台无关的。在我的Linux计算机上，sizeof(long)返回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来评估误差并没有改变结果太多。

- francis

---

感谢分析！ - alias