平方根倒数速算法

平方根倒数速算法（英語：Fast Inverse Square Root，亦常以“Fast InvSqrt()”或其使用的十六进制常数0x5f3759df代称）是用于快速计算${\displaystyle \textstyle x^{-1/2}}$（即${\displaystyle \textstyle x}$的平方根的倒数，在此${\displaystyle \textstyle x}$需取符合IEEE 754标准格式的32位浮点数）的一种算法。这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费，而在计算机图形学领域，若要求取照明和投影的波动角度与反射效果，就常需计算平方根倒数。

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

此算法最早被认为是由约翰·卡马克于90年代前期在SGI Indigo（英语：SGI Indigo）的开发中使用，后来则于1999年在《雷神之锤III竞技场》的源代码中应用，但直到2002－2003年间才在Usenet一类的公共论坛上出现^[1]。后来的调查显示，该算法在这之前就于计算机图形学的硬件与软件领域有所应用，如SGI和3dfx就曾在产品中应用此算法，但至今为止仍未能确切知晓算法中所使用的特殊常数的起源。

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

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

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

而${\displaystyle {\boldsymbol {\hat {v}}}={\boldsymbol {v}}/\|{\boldsymbol {v}}\|}$求得的就是标准化的向量，若以${\displaystyle {\boldsymbol {x}}}$代表${\displaystyle v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}$，则有${\displaystyle {\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（邪恶的浮点数位运算黑科技）
	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年代初（也即该算法发明的大概时间），软件开发时通用的平方根演算法（英语：Methods of computing square roots）多是从查找表中取得近似值^{[文 3]}，而这段代码取近似值耗时比之更短，达到精确度要求的速度也比通常使用的浮点除法计算法快四倍^{[文 4]}，虽然此算法会损失一些精度，但性能上的巨大优势已足以补偿损失的精度^{[文 5]}。由代码中对原数据的变量类型声明为float可看出，这一算法是针对IEEE 754标准格式的32位浮点数设计的，不过据Chris Lomont和后来的Charles McEniry的研究看，这一算法亦可套用于其他类型的浮点数上。

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

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

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

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

${\displaystyle m_{x}={\frac {M_{x}}{L}}}$

${\displaystyle e_{x}=E_{x}-B}$，其中${\displaystyle B=2^{b-1}-1}$

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

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

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

${\displaystyle \log _{2}{(y)}=-{\frac {1}{2}}\log _{2}{(x)}}$

${\displaystyle \log _{2}(1+m_{y})+e_{y}=-{\frac {1}{2}}\log _{2}{(1+m_{x})}-{\frac {1}{2}}e_{x}}$

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

${\displaystyle m_{y}+\sigma +e_{y}=-{\frac {1}{2}}m_{x}-{\frac {1}{2}}\sigma -{\frac {1}{2}}e_{x}}$

参照首段等式代入${\displaystyle M_{x}}$，${\displaystyle E_{x}}$，${\displaystyle B}$与${\displaystyle L}$，有：

${\displaystyle M_{y}+(E_{y}-B)L=-{\frac {3}{2}}\sigma {L}-{\frac {1}{2}}M_{x}-{\frac {1}{2}}(E_{x}-B)L}$

移项整理得：

${\displaystyle E_{y}L+M_{y}={\frac {3}{2}}(B-\sigma )L-{\frac {1}{2}}(E_{x}L+M_{x})}$

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

${\displaystyle I_{y}=E_{y}L+M_{y}=R-{\frac {1}{2}}(E_{x}L+M_{x})=R-{\frac {1}{2}}I_{x}}$，其中${\displaystyle R={\frac {3}{2}}(B-\sigma )L}$.

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

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

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

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

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

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

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

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

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

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

不仅该算法的原作者不明，人们也仍无法確定当初选择这个“魔术数字”的方法。Chris Lomont曾做了个研究：他推算出了一个函数以討論此速算法的誤差，並找出了使誤差最小的最佳R值0x5f37642f（与代码中使用的0x5f3759df相当接近），但以之代入算法计算并进行一次牛顿迭代后，所得近似值之精度仍略低於代入0x5f3759df的结果^{[文 17]}；因此Lomont将目标改为尋找在进行1－2次牛顿迭代后能得到最大精度的R值，在暴力搜尋後得出最优R值为0x5f375a86，以此值代入算法并进行牛顿迭代，所得的结果都比代入原始值（0x5f3759df）更精确^{[文 17]}，于是他的论文最后提到“如果可能我想询问原作者，此速算法是以数学推导还是以反复试错的方式求得？”^{[文 18]}。在论文中，Lomont亦指出，64位的IEEE754浮点数（即双精度类型）所对应的魔术数字是0x5fe6ec85e7de30da，但后来的研究表明，代入0x5fe6eb50c7aa19f9的结果精确度更高（McEniry得出的结果则是0x5FE6EB50C7B537AA，精度介于两者之间，Matthew Robertson给出的精度更高的值是0x5FE6EB50C7B537A9^[6]）。在Charles McEniry的论文中，他使用了一种类似Lomont但更复杂的方法来优化R值：他最开始使用穷举搜索，所得结果与Lomont相同^{[文 19]}；而后他尝试用調日法寻找最优值，所得结果恰是代码中所使用的魔术数字0x5f3759df，因此，McEniry认为，这一常数最初或许便是以“在可容忍误差范围内使用調日法”的方式求得^{[文 20]}。

• Kushner, David. The wizardry of Id. IEEE Spectrum. Aug 2002, 39 (8): 42–47. doi:10.1109/MSPEC.2002.1021943. 
• 顾森（Matrix67）. 旧闻一则：神秘的0x5f3759df 不可思议的Quake III源码. 2007-11-24 [2015-08-01]. （原始内容存档于2015-08-04）. 

• Origin of Quake3's Fast InvSqrt() （页面存档备份，存于互联网档案馆）, Beyond3D.com
• Quake III Arena source code （页面存档备份，存于互联网档案馆）, id Software
• Margolin, Tomer. Magical Square Root Implementation In Quake III. CodeMaestro. The Coding Experience. 2005-08-27. （原始内容存档于2012-04-14）.