0x5f3759df 是怎么算出来的?

(注:本文包含数学公式,使用 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

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

This entry was posted in Depth, Discussion, 中文 and tagged , , . Bookmark the permalink. Post a comment or leave a trackback: Trackback URL.

2 Comments

  1. Ahaxzh
    Posted December 14, 2011 at 3:56 pm | Permalink

    :)

    哇~发现您的网站支持IPv6访问呀!

Post a Comment

Your email is never published nor shared. Required fields are marked *

*
*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>