平方根倒数速算法

来源：互联网发布：数据备份厂商编辑：程序博客网时间：2024/04/27 20:54

维基百科，自由的百科全书

游戏实现光照和反射效果时以平方根倒数速算法计算波动角度，此图以第一人称射击游戏OpenArena为例。

平方根倒数速算法（英语：Fast Inverse Square Root，亦常以“Fast InvSqrt()”或其使用的十六进制常数0x5f3759df代称）是用于快速计算 $\scriptstyle x^{-\tfrac{1}{2}}$ （即 $\scriptstyle x$ 的平方根的倒数，在此 $\scriptstyle x$ 需取符合IEEE 754标准格式的32位浮点数）的一种算法。此算法最早可能是于90年代前期由SGI所发明，后来则于1999年在《雷神之锤III竞技场》的源代码中应用，但直到2002－2003年间才在Usenet一类的公共论坛上出现^[1]。这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费，而在计算机图形学领域来说求取照明和投影的波动角度与反射效果时就常需要计算平方根倒数。

此算法首先接收一个32位带符浮点数，然后将之作为一个32位整数看待，以将其向右进行一次逻辑移位的方式将之取半，并用十六进制“魔术数字”0x5f3759df减之，如此即可得对输入的浮点数的平方根倒数的首次近似值；而后重新将其作为浮点数，以牛顿法反复迭代以求出更精确的近似值，直至求出符合精确度要求的近似值。以此算法计算浮点数的平方根倒数的同一精度的近似值比直接使用浮点数除法要快四倍。

此算法最早被认为是由约翰·卡马克所发明，但后来的调查显示，该算法在这之前就于计算机图形学的硬件与软件领域有所应用，如SGI和3dfx就曾在产品中应用此算法。而就现在所知，此算法最早由Gary Tarolli在SGI Indigo的开发中使用。虽说在随后的相关研究中也提出了一些可能的来源，但至今为止仍未能确切知晓此常数的起源。

[编辑]算法的切入点

法线常在光影效果实现计算时使用，而这就涉及到矢量范数的计算。图中所标识的就是与一个面所垂直的一些矢量的集合。

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

要将一个矢量标准化，就必须计算其欧几里得范数以求得矢量长度，而这时就需对矢量的各分量的平方和求平方根；而当求取到其长度并以之除该矢量的每个分量后，所得的新矢量就是与原矢量同向的单位矢量，若以公式表示：

$\|\boldsymbol{v}\| = \sqrt{v_1^2+v_2^2+v_3^2}$ 可求得矢量v的欧几里得范数，此算法正类如对欧几里得空间的两点求取其欧几里得距离，

而 $\boldsymbol{\hat{v}} = \boldsymbol{v} / \|\boldsymbol{v}\|$ 求得的就是标准化的矢量，若以 $\boldsymbol{x}$ 代表 $v_1^2+v_2^2+v_3^2$ ，则有 $\boldsymbol{\hat{v}} = \boldsymbol{v} / \sqrt{x}$ ，

可见标准化矢量时需要用到对矢量分量的平方根倒数计算，所以对平方根倒数计算算法的优化对计算正规化矢量也大有裨益。

为了加速图像处理单元计算，《雷神之锤III竞技场》使用了平方根倒数速算法，而后来采用现场可编程逻辑门阵列的顶点着色器也应用了此算法^{[文 2]}。

[编辑]代码概览

下列代码是平方根倒数速算法在《雷神之锤III竞技场》源代码中的应用实例。示例剥离了C语言预处理器的指令，但附上了原有的注释^[2]：

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（对浮点数的邪恶位级hack）        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（第二次迭代，可以删除）         return y;}

为计算平方根倒数的值，软件首先要先确定一个近似值，而后则使用某些数值方法不断计算修改近似值，直至达到可接受的精度。在1990年代初（也即该算法发明的大概时间）软件里通用的平方根计算方法多是从查找表中取得近似值^{[文 3]}，而这段代码取近似值耗时比之更短，达到精确度要求的速度也比通常使用的浮点除法计算法快四倍^{[文 4]}，虽然此算法会损失一些精度，但性能上的巨大优势已足以补偿损失的精度^{[文 5]}。由代码中对原数据的变量类型声明为float可看出，这一算法是针对的是IEEE 754标准格式的32位浮点数设计的，不过据Chris Lomont和后来的Charles McEniry的研究来看，这一算法也可套用于其他类型的浮点数上。

平方根倒数速算法在速度上的优势源自将浮点数转化为长整型^{[注 1]}以作整数看待，并用特定常数0x5f3759df与之相减。然而对于代码阅读者来说，他们却难以立即领悟出使用这一常数的目的，因此这一常数和其它在代码中出现的难以理解的常数同被称为“魔术数字”^[1]^{[文 7]}^{[文 8]}^{[文 9]}。如此将浮点数当作整数先位移后减法所得的浮点数结果即是对输入数字的平方根倒数的粗略估计值，而后再进行一次牛顿迭代法以使之更精确后代码即执行完毕。由于算法所生成的用于输入牛顿法的首次近似值已经相当精确，所以此算法所计算出的近似值的精度已经可以接受，而若使用与《雷神之锤III竞技场》同为1999年发布的Pentium III中的SSE指令rsqrtss计算，则计算平方根倒数的收敛速度更慢，精度也更低^[4]。

[编辑]将浮点数转化为整数

参见：浮点数

要理解这段代码，首先需了解浮点数的存储格式。一个浮点数以32个二进制位表示一个有理数，而这32位由其意义分为三段：首先首位为符号位，如若是0则为正数，反之为负数；接下来的8位表示经过偏移处理（这是为了使之能表示-127－128）后的指数；最后23位表示的则是有效数字中除最高位以外的其余数字。将上述结构表示成公式即为 $\scriptstyle x=(-1)^{\mathrm{Si}}\cdot(1+m)\cdot 2^{(E-B)}$ ，其中 $\scriptstyle m$ 表示有效数字的尾数（此处 $\scriptstyle 0 \le m < 1$ ，偏移量 $\scriptstyle B=127$ ^{[文 10]}，而指数的值 $\scriptstyle E-B$ 决定了有效数字（在Lomont和McEniry的论文中称为“尾数”（mantissa））代表的是小数还是整数^{[文 11]}。以上图为例，将描述带入有 $\scriptstyle m=1\times 2^{-2}=0.250$ ），且 $\scriptstyle E-B=124-127=-3$ ，则可得其表示的浮点数为 $\scriptstyle x=(1+0.250)\cdot 2^{-3}=0.15625$ 。

符号位 01111111=12700000010=200000001=100000000=011111111=−111111110=−210000001=−12710000000=−128

8位二进制整数补码示例

如上所述，一个有符号正整数在二进制补码系统中的表示中首位为0，而后面的各位则用于表示其数值。将浮点数取别名存储为整数时，该整数的数值即为 $\scriptstyle I=E\times 2^{23}+M$ ，其中E表示指数，M表示有效数字；若以上图为例，图中样例若作为浮点数看待有 $\scriptstyle E=124$ ， $M=1\cdot 2^{21}$ ，则易知其转化而得的整数型号数值为 $I=124\times 2^{23} + 2^{21}$ 。由于平方根倒数函数仅能处理正数，因此浮点数的符号位（即如上的Si）必为0，而这就保证了转换所得的有符号整数也必为正数。以上转换就为后面的计算带来了可行性，之后的第一步操作（逻辑右移一位）即是使该数的长整形式被2所除^{[文 12]}。

[编辑]“魔术数字”

S(ign,符号)E(xponent,指数)M(antissa,尾数)1 位b位(n-1-b)位n位^{[文 13]}

对猜测平方根倒数速算法的最初构想来说，计算首次近似值所中使用的常数0x5f3759df也是重要的线索。为确定程序员最初选此常数以近似求取平方根倒数的方法，Charles McEniry首先检验了在代码中选择任意常数R所求取出的首次近似值的精度。回想上一节关于整数和浮点数表示的比较：对于同样的32位二进制数码，若为浮点数表示时实际数值为 $\scriptstyle x=(1+m_x)2^{e_x}$ ，而若为整数表示时实际数值则为 $\scriptstyle I_x=E_xL+M_x$ ^{[注 2]}，其中 $\scriptstyle L=2^{n-1-b}$ 。以下等式引入了一些由指数和有效数字导出的新元素：

$m_x=\frac{M_x}{L}$

$e_x=E_x-B$ ，其中 $B=2^{b-1}-1$

再继续看McEniry 2007里的进一步说明：

$y=\frac{1}{\sqrt{x}}$

对等式的两边取二进制对数（ $\textstyle \log_2$ ，即函数 $\textstyle f(n)=2^n$ 的反函数），有

$\log_2{(y)}=-\frac{1}{2}\log_2{(x)}$

$\log_2(1+m_y)+e_y=-\frac{1}{2}\log_2{(1+m_x)}-\frac{1}{2}e_x$

以如上方法就能将浮点数x和y的相关指数消去以将乘方运算化为加法运算。而由于 $\scriptstyle \log_2{(x)}$ 与 $\scriptstyle \log_2{(x^{-1/2})}$ 线性相关，因此在 $\scriptstyle x$ 与 $\scriptstyle y_0$ （即输入值与首次近似值）间就可以线性组合的方式创建方程^{[文 10]}。在此McEniry再度引入新数 $\sigma$ 描述 $\scriptstyle \log_2{(1+x)}$ 与近似值R间的误差^{[注 3]}：由于 $\scriptstyle 0 \le x < 1$ ，有 $\scriptstyle \log_2{(1+x)}\approx {x}$ ，则在此可定义 $\sigma$ 与x的关系为 $\scriptstyle \log_2{(1+x)}\cong x+\sigma$ ，这一定义就能提供二进制对数的首次精度值（此处 $0\le\sigma\le\tfrac{1}{3}$ ；当R为0x5f3759df时，有 $\scriptstyle \sigma=0.0450461875791687011756$ ^{[文 13]}）。由此将 $\scriptstyle \log_2{(1+x)}= x+\sigma$ 代入上式，有：

$m_y+\sigma+e_y=-\frac{1}{2}m_x-\frac{1}{2}\sigma-\frac{1}{2}e_x$

参照首段等式代入 $M_x$ ， $E_x$ ， $B$ 与 $L$ ，有：

$M_y+(E_y-B)L=-\frac{3}{2}\sigma{L}-\frac{1}{2}M_x-\frac{1}{2}(E_x-B)L$

移项整理得：

$E_yL+M_y=\frac{3}{2}(B-\sigma)L-\frac{1}{2}(E_xL+M_x)$

如上所述，对于以浮点规格存储的正浮点数x,若将其作为长整型表示则示值为 $\textstyle I_x=E_xL+M_x$ ，由此即可根据x的整数表示导出y（在此 $\textstyle y=\frac{1}{\sqrt{x}}$ ，亦即x的平方根倒数的首次近似值）的整数表示值，也即：

$I_y=E_yL+M_y=R-\frac{1}{2}(E_xL+M_x)=R-\frac{1}{2}I_x$ ，其中 $R=\frac{3}{2}(B-\sigma)L$ .

最后导出的等式 $\scriptstyle I_y=R-\frac{1}{2}I_x$ 即与上节代码中i = 0x5f3759df - (i>>1);一行相契合，由此可见平方根倒数速算法中对浮点数进行一次移位操作与整数减法就可以可靠地输出一个浮点数的对应近似值^{[文 13]}。到此为止McEniry的证明只显示出了可用常数R为工具以近似求取浮点数的平方根倒数，但仍未能确定代码中的R值的选取方法。

关于作一次移位与减法操作以使浮点数的指数被-2除的原理，Chris Lomont的论文中亦有个相对简单的解释：以 $\scriptstyle 10000=10^4$ 为例，将其指数除-2可得 $\scriptstyle 10000^{-1/2}=10^{-2}=1/100$ ；而由于浮点表示的指数有进行过偏移处理，所以指数的真实值e应为 $\scriptstyle e=E-127$ ，因此可知除法操作的实际结果为 $\scriptstyle -e/2+127$ ^{[文 14]}，这时用R（在此即为“魔术数字”0x5f3759df）减之即可使指数的最低有效数位转入有效数字域，之后重新转换为浮点数时，就能得到一个相当接近所输入的浮点数的平方根倒数的近似值。在这里对常数R的选取亦有所讲究，选取一个好的R值可以减少对指数进行除法与对有效数字域进行移位时可能产生的错误。基于这一标准，0xbe即是最合适的R值，而0xbe右移一位即可得到0x5f，这恰是魔术数字R的第一个字节^{[文 15]}。

[编辑]精确度

使用启发式平方根倒数速算法与使用C语言标准库libstdc的函数所计算出的平方根倒数的差值一览，注意这里使用的是双对数坐标系。

如上所述，平方根倒数速算法所得的近似值惊人的精确，右图亦展示了以上述代码计算（以平方根倒数速算法计算后再进行一次牛顿法迭代）所得近似值的误差：当输入0.01时，以C语言标准库函数计算可得10.0，而InvSqrt()得值为9.9825822，其间误差为0.017479，相对误差则为0.175%，且当输入更大的数值时，绝对误差不断下降，相对误差也一直控制在一定的范围之内。

[编辑]牛顿法

主条目：牛顿法

在进行了如上的整数操作之后，示例程序再度将被转为长整型的浮点数回转为浮点数（对应x = *(float*)&i;）并对其进行一次浮点运算操作（对应x = x*(1.5f - xhalf*x*x);），这里的浮点运算操作就是对其进行一次牛顿法迭代，若以此例说明：

$y=\frac{1}{\sqrt{x}}$ 所求的是y的平方根倒数，以之构造以y为自变量的函数，有 $f(y)=\frac{1}{y^2}-x=0$ ，

将其代入牛顿法的通用公式 $y_{n+1} = y_{n} - \frac{f(y_n)}{f'(y_n)}$ （其中 $\, y_n$ 为首次近似值），

有 $y_{n+1} = \frac{y_{n}(3-xy_n^2)}{2}$ ，其中 $f(y)=\frac{1}{y^2}-x$ ， $f'(y)=\frac{-2}{y^3}$ 。

整理有 $\, y_{n+1} = \frac{y_{n}(3-xy_n^2)}{2} = y_{n}(1.5-\frac{xy_n^2}{2})$ ，对应的代码即为x = x*(1.5f - xhalf*x*x);。

在以上一节的整数操作产生首次近似值后，程序会将首次近似值作为参数送入函数最后两句进行精化处理，代码中的两次迭代（以一次迭代的输出（对应公式中的 $y_{n+1}$ ）作为二次迭代的输入）正是为了进一步提高结果的精度^{[文 16]}，但由于雷神之锤III引擎的图形计算中并不需要太高的精度，所以代码中只进行了一次迭代，二次迭代的代码则被注释^{[文 9]}。

[编辑]历史与考究

id Software的创始人约翰·卡马克。虽然这段代码并非他所作，但他常被认为与这段代码相关。

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

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

[编辑]注释

^ 由于现代计算机系统对长整型的定义有所差异，使用长整型会降低此段代码的可移植性。具体来说，由此段浮点转换为长整型的定义可知，如若这段代码正常运行，所在系统的长整型长度应为4字节（32位），否则重新转为浮点数时可能会变成负数；而由于C99标准的广泛应用，在现今多数64位计算机系统（除使用LLP64数据模型的Windows外）中，长整型的长度都是8字节^{[文 6]}^[3]。
^ 此处“浮点数”所指为标准化浮点数，也即有效数字部分必须满足 $\scriptstyle 0 \le m_x <1)$ ，可参见David Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys. 1991.March, 23(1): 5–48. doi:10.1145/103162.103163.
^ Lomont 2003确定R的方式则有所不同，他先将R分解为 $\scriptstyle R_1$ 与 $\scriptstyle R_2$ ，其中 $\scriptstyle R_1$ 与 $\scriptstyle R_2$ 分别代表R的有效数字域和指数域^{[文 7]}。

[编辑]参考

[编辑]脚注

^ ^1.0 ^1.1 ^1.2 ^1.3 Sommefeldt, Rys. Origin of Quake3's Fast InvSqrt(). Beyond3D. 2006-11-29 [2009-02-12] （英文）.
^ quake3-1.32b/code/game/q_math.c. Quake III Arena. id Software. [2010-11-15].
^ Meyers, Randy. The New C: Integers in C99, Part 1. drdobbs.com. 2000-12-01 [2010-09-04].
^ Ruskin, Elan. Timing square root. Some Assembly Required. 2009-10-16 [2010-09-13].
^ Sommefeldt, Rys. Origin of Quake3's Fast InvSqrt() - Part Two. Beyond3D. 2006-12-19 [2008-04-19].