c - Quake 反平方根 : accuracy

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

paper 声明对于所有正正态 float ，相对误差最多为 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

因此，这个特定的输入似乎违反了该论文中的声明。

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

Answer 1

您使用了错误的魔数(Magic Number)。

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