卡马克快速平方根倒数算法

逛知乎的时候看到这么个问题:

其中有个答案提到了一种计算浮点数的平方根倒数的快速算法,其实我以前就曾经思考过大数平方根的计算方法,然后…就没有然后了。

神奇的地方在于这个快速算法中有个诡异的常数,据说至今没有人知道它最早是怎么来的。

代码示例

下面这段代码来源于上面那个答案。

1
2
3
4
5
6
7
8
9
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;
}

这个算法说到头只有这么几行而已,而恰恰神奇的就是这么几行代码可以非常快速地实现平方根倒数。

整个算法的关键就是第5行的那个诡异常数以及第7行的牛顿迭代了。

算法原理

要理解这个算法,首先要明白牛顿迭代,而牛顿迭代法的基础是泰勒级数

泰勒级数

定义:如果$f(x)$在点$x=x_0$处具有任意阶导数,则幂级数:

$$
\begin{align}
f(x)&=\sum^{\infty}_{n=0}\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n\\
&=f(x_0)+f’(x_0)(x-x_0)+\frac{f’’(x_0)}{2!}(x-x_0)^2+…+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+…
\end{align}$$

称为$f(x)$在$x=x_0$处的泰勒级数。

泰勒展开在数学上还有很多的应用…然而,高数学过去那么多年了,有啥用我早忘了…0.0

这里只要知道泰勒级数可以用来近似计算函数值就行了。

后面部分明显是越往后越小的,收敛性应该是要有条件来严格证明的…这个我也忘了。

牛顿迭代

牛顿迭代解线性方程,就是把线性方程$f(x)=0$线性化的一种近似方法。把$f(x)$在$x_0$的某个邻域内展开成泰勒级数,取其线性部分(前两项):

$$
f(x)=f(x_0)+f’(x_0)(x-x_0)
$$

令其为零,就可以作为原本$f(x)=0$方程的近似方程了。只要$f’(x_0)\neq0$,方程的解就可以写成:

$$
x_1=x_0-\frac{f(x_0)}{f’(x_0)}
$$

这样,就得到牛顿迭代法的迭代关系式了:

$$
x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)}
$$

理论上来说,是不是我任意取一个满足$f’(x_0)\neq0$的数,都可以把它作为第1个近似解,然后通过不断地进行牛顿迭代,直到精度达到要求即解出了比较靠谱的方程解?


回到目标算法上面来,假设已知的数是$a$,要求它的平方根的倒数,就是解方程:

$$
(\frac{1}{x})^2=a
$$

即解方程:

$$
f(x)=\frac{1}{x^2}-a=0
$$

那么:

$$
f’(x)=-\frac{2}{x^3}
$$

代回到牛顿迭代式中去:

$$
\begin{align}
x_{n+1}&=x_n-\frac{f(x_n)}{f’(x_n)}\\
&=x_n-\frac{\frac{1}{x^2_n}-a}{-\frac{2}{x^3_n}}\\
&=x_n+\frac{x_n-ax^3_n}{2}\\
&=\frac{3}{2}x_n-\frac{a}{2}x^3_n
\end{align}
$$

这个就是上面程序第7行的那个牛顿迭代式了。上面那个程序段只迭代了一次,其实精度是不够的,一般迭代2次,跟标准值的误差就不大了。

当然,迭代次数越多,精度越高,相应地耗时也是越长。

诡异的常数

最后,还有个问题没解决……0x5f3759df这个常数是什么鬼啊?

上面已经讨论过了,通过牛顿迭代法,我们可以在第n个近似解的基础上推出第n+1个精度更高的近似解。这里唯一的漏洞就是:我们必须得到第1个近似解!这样才能够往下递推!

而上面第5行的这段代码:

1
i = 0x5f3759df - (i >> 1);

就是用来计算第1个近似解的。


好吧,还是不明白这个数是怎么来的…

据史料记载:

浮点数的平方根倒数常用于计算正规化矢量。3D图形程序需要使用正规化矢量来实现光照和投影效果,因此每秒都需做上百万次平方根倒数运算,而在处理坐标转换与光源的专用硬件设备出现前,这些计算都由软件完成,计算速度亦相当之慢;在1990年代这段代码开发出来之时,多数浮点数操作的速度更是远远滞后于整数操作,因而针对正规化矢量算法的优化就显得尤为重要。

1999年的《雷神之锤III竞技场》就借助了这个神奇的算法。

《雷神之锤III》的代码直到QuakeCon 2005才正式放出,但早在2002年(或2003年)时平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了。最初人们猜测是卡马克写下了这段代码,但他在询问邮件的回复中否定了这个观点,并猜测可能是先前曾帮id Software优化雷神之锤的资深汇编程序员Terje Mathisen写下了这段代码;而在Mathisen的邮件里他表示在1990年代初他只曾作过类似的实现,确切来说这段代码亦非他所作。现在所知的最早实现是由Gary Tarilli在SGI Indigo中实现的,但他亦坦承他仅对常数R的取值做了一定的改进,实际上他也不是作者。Rys Sommefeldt则在向以发明MATLAB而闻名的Cleve Moler查证后认为原始的算法是Ardent Computer公司的Greg Walsh所发明,但他也没有任何决定性的证据能证明这一点。

目前不仅该算法的原作者不明,人们也仍无法明确当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验:他编写了一个函数,以在一个范围内遍历选取R值的方式将逼近误差降到最小,以此方法他计算出了线性近似的最优R值0x5f37642f(与代码中使用的0x5f3759df相当接近),但以之代入算法计算并进行一次牛顿迭代后,所得近似值与代入0x5f3759df的结果相比精度却仍略微更低;而后Lomont将目标改为遍历选取在进行1-2次牛顿迭代后能得到最大精度的R值,并由此算出最优R值为0x5f375a86,以此值代入算法并进行牛顿迭代后所得的结果都比代入原始值(0x5f3759df)更精确,于是他的论文最后以“原始常数是以数学推导还是以反复试错的方式求得”的问题作结。在论文中Lomont亦指出64位的IEEE754浮点数(即双精度类型)所对应的魔术数字是0x5fe6ec85e7de30da,但后来的研究表明代入0x5fe6eb50c7aa19f9的结果精确度更高(McEniry得出的结果则是0x5FE6EB50C7B537AA,精度介于两者之间)。在Charles McEniry的论文中,他使用了一种类似Lomont但更复杂的方法来优化R值:他最开始使用穷举搜索法,所得结果与Lomont相同;而后他尝试用带权二分法寻找最优值,所得结果恰是代码中所使用的魔术数字0x5f3759df,因此McEniry确信这一常数或许最初便是以“在可容忍误差范围内使用二分法”的方式求得。


先复习一下计算机中浮点数的表示法

然后回到上面讨论过的第1个近似解,第1个近似解越接近正确答案,则显然牛顿迭代法的效果会越好。

考虑已知的浮点数$a$在计算机中的实际表示为:

$$(1+F)*2^E$$

那么它的平方根倒数为:

$$(1+F)^{-\frac{1}{2}}*2^{-\frac{E}{2}}$$

猜测一下:第5行代码中,将i右移1位,再用常数减掉,可以看成是将指数除2之后再取反。至于尾数…早都变得不知道哪去了。


所以,还是没明白这个常数是怎么来的,我只能猜测到这里了。

这里把Chris Lomont写的那篇论文也一并贴出来:FAST INVERSE SQUARE ROOT

0%