数学之美：平方根倒数速算法中的神奇数字 0x5f3759df

http://blog.jobbole.com/105295/?ref=myread

本文由 伯乐在线 - JLee 翻译，黄利民 校稿
英文出处：Christian Plesner Hansen。

【伯乐在线导读】：在伯乐小组的这个讨论帖《你见过或写过的最复杂的 C 语言程序是？》中，就提到了「平方根倒数速算法」中的神奇数字 0x5f3759df。Christian 在本文中探讨了该算法诸多有趣的地方，解释了它的原理，对指数在 -1 到 1 的情况进行了拓展，还有一些数学相关的新思路。

0x5f3759df

这篇文章讲述一个神奇的常数 0x5f3759df 和一个非常简洁的算法，平方根倒数速算法，也是这个常数的来源。

来看看平方根倒数速算法：

C

1

2

3

4

5

6

7

8

float FastInvSqrt(float x) {

  float xhalf = 0.5f * x;

  int i = *(int*)&x;         // evil floating point bit level hacking

  i = 0x5f3759df - (i >> 1);  // what the fuck?

  x = *(float*)&i;

  x = x*(1.5f-(xhalf*x*x));

  return x;

}

这段代码可以快速计算

的高精度近似值。

现如今这个函数很出名，而它的流行是从出现在 2005 年的 Quake III Arena 的源码中开始的，起初归功于 John Carmack，但后来发现早在 Quake 之前，80 年代中期，Ardent Computer 公司就已使用，比 SGI 和 3dfx 还早，已能追溯到的最初作者是 Greg Walsh。上面具体的代码是 Quake 源码的改编版（Quake 源码是所有讨论的源头）。

这篇文章探讨了这个算法诸多有趣的地方，解释了它的原理，对指数在 -1 到 1 的情况进行了拓展，还有一些数学相关的新思路。

（本文确实包含许多数学知识，你可以把这些方程看做注释，无需阅读即可掌握本文要点，但如果你想要完整的阅读本文，验证我的观点的正确性，则需要关注这些方程。）

为什么？

为何需要计算倒平方根，甚至于值得实现一个怪异的算法来加速计算？因为它一直是 3D 编程计算中的一部分。在 3D 图形中，你使用平面法线，长度为 1 的三坐标向量，来表示光线和反射。你会使用很多平面法线，计算它们时需要对向量进行标准化。如何标准化一个向量呢？每个坐标分别除以向量的长度，因此，每个坐标需乘上

计算  相对来说开销很小，计算平方根并做除法开销很大。进入FastInvSqrt。

如何做到的？

这个函数究竟做了什么计算出了结果？它有 4 个主要步骤。首先它把浮点类型输入的二进制表示重定义为一个整形数字的二进制表示。

C

1

inti=*(int*)&x;        // evil floating point bit level hack

运用得到的结果进行整形算术计算，得到我们所求值的近似值:

C

1

i = 0x5f3759df - (i >> 1);  // what the fuck?

尽管结果不是近似值本身，但如果你把这个整形数字的二进制表示重定义为一个浮点数，那么就会得到近似值。所以代码做了与步骤一相反的变换回到了浮点数：

C

1

x=*(float*)&i;

最终进行一次牛顿法迭代来提高精度。

C

1

x = x*(1.5f-(xhalf*x*x));

这便得到 x 倒平方根的非常好的近似值。最后一步牛顿法相对比较直观，我便不多花时间讲解。关键步骤是第二步：对原始浮点数转化来的整形数进行算术运算，得到一个有意义的结果，下面重点关注这一步。

怎么回事？

这部分解释第二步背后的数学原理（下面推导的第一部分，至常数值的计算，最早由 McEniry 发现的）。

在进入这一有趣的部分之前，我先快速解释下标准浮点数编码，我只讲解我要用的部分，完整的背景知识请参见维基百科。浮点数由三部分组成：符号，指数和尾数。这是单精度（32 位）浮点数二进制表示：

C

1

seeeeeeeemmmmmmmmmmmmmmmmmmmmmmm

最高位是符号，接下来的 8 位是指数，底部 23 位是尾数。由于将要计算的平方根只对正数定义，所以从现在起假设符号位是 0 。

当把一个浮点数看作一堆 0 、1 时，指数和底数便只是普通的正整数，没什么特别的。把它们记为 E 和 M（会经常使用）。另一方面，当把这些 0、1 解释为浮点值时，尾数则是介于 0 到 1 之间的值，全 0 意味着 0，全 1 是非常接近但略微小于 1 的数。并非将指数当做 8 位无符号整形使用，而是减去一个偏量 B，使之成为介于 -127 到 128 之间的有符号整形。把这些值的浮点解释记为 e 和 m。总之，我按照 McEnry 把整形解释相关的值用大写字母表示，把浮点解释相关的值用小写字母表示。

这两种解释之间的转换很直接：

对 32 位浮点数来说， L 是 223，B 是 127。给定 e 和 m，浮点数的值可计算得到：

相应地整形解释的数是

现在解释这个算法的所有基础几乎都有了，因此可以开始了，遗漏的部分会在解释过程中引入。给定输入 x，想要计算的是它的倒平方根，或者

我们先对等式两边取以 2 为底的对数，至于原因很快就会知道：

既然我们进行计算的值都是浮点数，则可以把 x 和 y 用各自的浮点组成替换：

对数比较麻烦，幸运的是可以很轻松的摆脱它们，但在这之前，需暂停一下来解释其原理。

在方程的两边都含有这样的项：

其中 v 介于 0 到 1 之间，也仅仅当 v 在 0 到 1 之间时，这个函数和一条直线很接近：

或者方程形式：

其中 σ 是一个可选常数，尽管不能完美匹配但可以调节 σ 使得二者非常接近。通过它我们可以把上述包含对数的方程近似为一个线性方程，计算起来更轻松：

到达关键的部分了！这时，用整形解释下的指数和尾数替代浮点解释下的表示进行计算就很方便了：

通过少许步骤进行移项合并，就会得到很熟悉的形式（细节很乏味，可以跳过）：

走完最后一步之后，有趣的事情发生了：在这些杂乱的项中，方程两边都有整形解释的值：

换句话说，y 的整形解释等于某个常数减去 x 的整形解释的一半，或者用 C 代码表示：

C

1

i = K - (i >> 1);

K 看起来很熟悉对吧？

剩下的工作就是找到这个常数。我们已经知道 B 和 L 的值，σ 的值还不知道。记住，σ 是对对数函数近似的调节因子，所以它的取值有些自由。我选取最初应用过的一个值，0.0450465，得到：

想知道这个值的 16 进制表示？0x5f3759df。（当然理应是它，因为选的特定的 σ 得到的），这个常数不是一个位模式，这点你可能通过它被写成 16 进制形式而推知，它只是一个普通计算的结果化为整形的形式。

但是正如 Knuth 所说，目前为止我们只是证明这个方法应该有效，还没检测。为了对这个近似的精确程度有直观的认识，用它的图线和精确的倒平方根的图线进行对比：

这里输入的值介于 0 到 100 之间。很精确对吗？理应如此。这不仅仅是神奇，正如上面介绍的其推导，它的计算也很奇怪，但都是清晰、有意义的操作——浮点和整形之间的位转换

等等，还有更多！

这个操作的推导告诉你的不只是一个常数的值，你会注意到这个推导不依赖于任何项的具体值，它们只是用于变换的常数，这意味着即使我们改变它们，推导也成立。

首先，这个计算不关心 L 和 B 具体是什么，它们通过浮点表示给出，这意味着对 64 位和 128 位浮点数可以使用同样的方法，所有要做的只是重新计算常数，这是唯一依赖它们的部分。

第二，它不关心 σ 的取值，使得对数函数和 x+σ 之间差异最小的 σ 可能不会得到最精确的近似，事实确实如此。这是浮点数的舍入和牛顿法联合造成的，σ 取值本身是一个有趣的课题，McEniry 和 Lomont 有过论述。

最后，不仅仅只能计算倒平方根，这里指数碰巧是 -1/2，但是对于 -1 到 1 之间的指数，推导依然成立。我们把指数记为 p（因为 e 已经使用），替代掉 -1/2 作同样的推导，得到：

来试试其他一些指数，首先 p=0.5，普通平方根计算：

或者代码形式，

C

1

i=0x1fbd1df5+(i>>1);

同样有效吗？当然：

这也许是近似计算平方根的一个很著名的方法，但 Google 和维基百科的一个粗略调查显示它并不是。

甚至对奇数次开方也适用，如立方根

相应地：

C

1

i = (int) (0x2a517d3c + (0.333f * i));

由于这个奇数因子，我们无法使用移位来代替乘法。同样近似很接近：

这时，你可能会注意到，当指数改变时，我们只需做非常简单的工作：只是通过改变一个线性因子来调整常数值，改变和输入的整形解释相乘的因子。这些操作开销都不大，因此在运行时进行是可行的，而不需要重新计算。如果我们事先对其他两个因子的乘积进行计算：

则不用提前知道指数，也可计算：

C

1

i=(1-p)*0x3f7a3bea+(p*i);

对项进行整理，还可以减少一步乘法计算：

C

1

i = 0x3f7a3bea + p * (i - 0x3f7a3bea);

这便是对 -1 到 1 之间任意指数进行快速指数计算的神奇部分，现在我们还差一部分，就是对牛顿近似进行一般化，使得快速指数计算函数对所有指数都适用，且和之前的倒平方根函数一样准确。这部分我没深入，就交给其他博客文章了（很可能是其他人的而不是我）。

上面表达式包含了一个新的神奇的常数，  0x3f7a3bea。即使它在某种程度上比最初的常数更神奇，最初的常数可由它演变而来，但它依赖于 σ 的主观选取，所以它也不具通用性。把这个常数记为 Cσ ，马上我们会更进一步探究它。

但首先，我们可以对 p=0 时的公式进行合理性检验。对于 p 取零时，结果应该始终为 1，因为 x0 等于 1，与 x 无关。确实，第二项消失了，因为它乘以了 0 ，我们得到简洁的：

C

1

i=0x3f7a3bea;

确实是常数，当解释为浮点值时等于 0.977477，近似为 1，所以合理性检验成立。这也告诉我们一些信息： 当 Cσ  转化为浮点解释时，是一个有意义的值，为 1 或非常接近 1。

这很有趣，更进一步，Cσ 的整形解释是：

这很像但不是一个浮点数的形式，唯一的问题是这里是减去而不是加上第二项，这个问题可以很轻易解决：

现在它看起来像一个浮点数的整形解释，为了具体来看，我们先分辨出指数和尾数，然后计算浮点数值，cσ, 这是指数：

这是尾数：

所以常数的浮点数值是：

确实如果你先用 σ 作除法，0.0450465 除以 2 得 0.02252325，用 1 减去它的 0.97747675 或者刚才的“近似为 1”。这让我们从另一种角度来看待 Cσ，作为一个浮点数的整形解释。用代码计算它：

C

1

2

3

floatsigma=0.0450465;

floatc_sigma=1-(0.5f*sigma);

intC_sigma=*(*int)&c_sigma;

注意，对于固定的 σ , 这些皆为常数，编译器应该能够优化整个计算过程。结果为 0x3f7a3beb，不是之前的 0x3f7a3bea ，但只有一字节的不同（最次要的一个字节），这对浮点数计算很正常。要得到最初的常数，本文的标题，只需把结果乘以 1.5。

我们已经足够深入了，至少我很满足了，没有什么可以继续探究的了。对于我来说，通过这个练习，主要学到的是整形和浮点型之间的位转换操作不是没有意义的，它很怪异但却很经济的数字操作，在计算中很有用。我认为它还有更多有待发现的用处。