平方根倒数速算法

✍ dations ◷ 2025-02-25 05:37:49 #平方根倒数速算法

平方根倒数速算法（英语：Fast Inverse Square Root，亦常以“Fast InvSqrt()”或其使用的十六进制常数0x5f3759df代称）是用于快速计算 ${displaystyle textstyle x^{-1/2}}$ ））代表的是小数还是整数。以上图为例，将描述带入有 ${displaystyle textstyle m=1times 2^{-2}=0.250}$ ${displaystyle textstyle m=1times 2^{-2}=0.250}$ ），且 ${displaystyle textstyle E-B=124-127=-3}$ ${displaystyle textstyle E-B=124-127=-3}$ ，则可得其表示的浮点数为 ${displaystyle textstyle x=(1+0.250)cdot 2^{-3}=0.15625}$ ${displaystyle textstyle x=(1+0.250)cdot 2^{-3}=0.15625}$ 。

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

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

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

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

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

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

移项整理得：

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

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

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

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

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

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

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

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