typeof.net

（注：本文包含数学公式，使用 SVG 图形格式，请不要在 RSS 阅读器里阅读本文）

基本上每个程序员都会对 Quake 那个神一样的常量 0x5f3759df 感到惊奇。在 Quake 里，有一段代码是计算光照时候使用的函数，计算 ，那段代码写的真叫经典：

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

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;  // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 ); // WTF?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

  return y
}

任何学过数值分析的人肯定都知道，最后两行干的事情是牛顿迭代法，用迭代逼近算出准确值。整个过程里，先把单精度浮点 number “硬”转换成整数，进行一个非常诡异的减法后捣腾回实数，然后进行牛顿迭代。那个诡异的步骤显然是算出 的一个近似值，而且从只迭代了一次看，猜测的非常好。因为是计算光照用，千分之一的误差是可以允许的，否则，多迭代几行就是。

这肯定会引来无数人遐想的。不过你看过下面一段内容后就知道其中就里了。

根据 IEEE-754，每一个正的二进制浮点数都可以表示成 。其中， 为其小数位，而 是其指数部分。对于单精度来说， 有 23 个二进制位，而 则有 8 位， 。

我们把它转换成“无符号整数形式” 。对比两者可以得到：

那好，为了计算 ，两边取对数，有：

为了去掉那个该死的对数，考虑到 在[0, 1] 上， ，在这里引入一个常数 σ 使得 。于是：

这样就可以解释那一行里那个著名的 0x5f3759df 实际上就是 。取 σ = 0 代入尝试，可得到 R = 1598029824 = 0x5f400000。这个数值已经比较接近了，然而还不准确，毕竟这里 σ = 0。

那么 σ 设置成多少合适呢？由于希望用牛顿法计算，自然需要让 和 差异尽量小。我使用的是最小二乘法，即计算积分：

积分 的结果是 σ 的二次函数，可算出其最小值时，σ 为 。利用这个 σ 计算出的 R = 0x5f34ff58。Quake 里使用的 σ = 0.0450461875791687011756，比用最小二乘法算出的 σ 小。另一个可以选用的 σ 是 在 [0, 1] 上的最大值的一半，为 ，算出的 R 为 1597488340 = 0x4f37bcb6。

当然了，用最小二乘法去“猜” σ 并不是非常好，最好的 σ 肯定是用大量实验选择的， 这个 σ 值并不是非常好的值，但是作为一次尝试，已经是很不错了。