0x5f3759df的数学原理 ----一种快速开方根求倒原理

0x5f3759df的数学原理 ----一种快速开方根求倒原理

转载于：

1.0x5f3759df的数学原理

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

3.0x5f3759df

Quake-III Arena (雷神之锤3)是90年代的经典游戏之一。

 

该系列的游戏不但画面和内容不错，而且即使计算机配置低，也能极其流畅地运行。这要归功于它3D引擎的开发者约翰-卡马克（John Carmack）。事实上早在90年代初DOS时代，只要能在PC上搞个小动画都能让人惊叹一番的时候，John Carmack就推出了石破天惊的Castle Wolfstein, 然后再接再励，doom, doomII, Quake...每次都把3-D技术推到极致。他的3D引擎代码资极度高效，几乎是在压榨PC机的每条运算指令。当初MS的Direct3D也得听取他的意见，修改了不少API。

最近，QUAKE的开发商ID SOFTWARE遵守GPL协议，公开了QUAKE-III的原代码，让世人有幸目睹Carmack传奇的3D引擎的原码。

我们知道，越底层的函数，调用越频繁。3D引擎归根到底还是数学运算。

那么找到最底层的数学运算函数（game/code/q_math.c），必然是精心编写的。里面有很多有趣的函数，很多都令人惊奇，估计我们几年时间都学不完。

在/code/game/q_math.c里发现了这样一段代码。

它的作用是将一个数开平方并取倒，经测试这段代码比(float)(1.0/sqrt(x))快4倍：

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 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

#ifndef Q3_VM
#ifdef __linux__
	assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
	return y;
}

函数返回1/sqrt(x)，这个函数在图像处理中比sqrt(x)更有用。

注意到这个函数只用了一次叠代！（其实就是根本没用叠代，直接运算）。编译，实验，这个函数不仅工作的很好，而且比标准的sqrt()函数快4倍！要知道，编译器自带的函数，可是经过严格仔细的汇编优化的啊！

这个简洁的函数，最核心，也是最让人费解的，就是标注了“what the fuck?”的一句：i = 0x5f3759df - ( i >> 1 );

再加上y = y * ( threehalfs - ( x2 * y * y ) );

两句话就完成了开方运算！而且注意到，核心那句是定点移位运算，速度极快！特别在很多没有乘法指令的RISC结构CPU上，这样做是极其高效的。

算法的原理其实不复杂,就是牛顿迭代法,用x-f(x)/f'(x)来不断的逼近f(x)=a的根。

简单来说比如求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入x-f(x)/f'(x)后有(x+a/x)/2，现在我们选a=5,选一个猜测值比如2，那么我们可以这么算

5/2 = 2.5;

(2.5+2)/2 = 2.25;

5/2.25 = xxx;

(2.25+xxx)/2 = xxxx

...

这样反复迭代下去，结果必定收敛于sqrt(5)，没错，一般的求平方根都是这么算的，但是卡马克(quake3作者)真正牛B的地方是他选择了一个神秘的常数0x5f3759df 来计算那个猜测值，就是我们加注释的那一行,那一行算出的值非常接近1/sqrt(n),这样我们只需要2次牛顿迭代就可以达到我们所需要的精度.好吧 如果这个还不算NB,接着看:

普渡大学的数学家Chris Lomont看了以后觉得有趣，决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。Lomont也是个牛人，在精心研究之后从理论上也推导出一个最佳猜测值，和卡马克的数字非常接近, 0x5f37642f。卡马克真牛，他是外星人吗？

传奇并没有在这里结束。Lomont计算出结果以后非常满意，于是拿自己计算出的起始值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确的求得平方根。结果是卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。

最后Lomont怒了，采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数字要好上那么一丁点的数字，虽然实际上这两个数字所产生的结果非常近似，这个暴力得出的数字是0x5f375a86。

Lomont为此写下一篇论文，"Fast Inverse Square Root"。

最后，给出最精简的1/sqrt()函数：

float InvSqrt(float x)
{
    float xhalf = 0.5f * x;
    int i = *(int *)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float *)&i;
    x = x * (1.5f - xhalf * x * x);
    return x;
}

为什么？

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

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

如何做到的？

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

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

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

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

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

x = *(float*)&i;

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

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

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

怎么回事？

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

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

s e e e e e e e e m m m m m m m m m m m m m m m m m m m m m m m

最高位是符号，接下来的 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 代码表示：

	
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，普通平方根计算：

或者代码形式，

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

同样有效吗？当然：

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

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

相应地：

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

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

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

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

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

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

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

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

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

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

i = 0x3f7a3bea;

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

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

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

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

这是尾数：

所以常数的浮点数值是：

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

float sigma = 0.0450465;
float c_sigma = 1 - (0.5f * sigma);
int C_sigma = *(*int)&c_sigma;

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

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