平方根倒数速算法

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

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

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

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

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

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年代初（也即该算法发明的大概时间），软件开发时通用的多是从查找表中取得近似值，而这段代码取近似值耗时比之更短，达到精确度要求的速度也比通常使用的浮点除法计算法快四倍，虽然此算法会损失一些精度，但性能上的巨大优势已足以补偿损失的精度。由代码中对原数据的变量类型声明为float可看出，这一算法是针对IEEE 754标准格式的32位浮点数设计的，不过据Chris Lomont和后来的Charles McEniry的研究看，这一算法亦可套用于其他类型的浮点数上。

平方根倒数速算法在速度上的优势源自将浮点数转化为长整型以作整数看待，并用特定常数0x5f3759df与之相减。然而对于代码阅读者来说，他们却难以立即领悟出使用这一常数的目的，因此和其它在代码中出现的难以理解的常数一样，这一常数亦被称为“魔術數字}-”。

未定义行为的问题
从 C99 和 C++98 标准开始，由于严格别名规则，形如 i = (long) &y; 的写法是未定义行为。想要正确做到类型双关，一种同时符合标准 C 和 C++ 的方式是使用 memcpy，但这不是唯一一种方法。

下面是更符合标准 C 和标准 C++，可移植性更高的代码：

#include
#include

float Q_rsqrt(float number)
{
float x = number;

// 第一步：将 x 当作整数，进行移位和减法操作，变成近似的结果
uint32_t i;
memcpy(&i, &x, sizeof i); // 编译器往往是会优化的，不必担心时间和空间开销
i = 0x5f3759dfu - (i >> 1);
memcpy(&x, &i, sizeof x);

// 第二步：对 x 进行牛顿迭代，使之更接近结果
const float xhalf = number * 0.5f;
x = x (1.5f - (xhalf x * x));
//x = x (1.5f - (xhalf x * x)); 第二次迭代可以去掉了

return x;
}

但下文仍然会用第一个版本的代码进行讲解。

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

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

“-{zh-hans:魔术数字;zh-hant:魔術數字}-”
对猜测平方根倒数速算法的最初构想来说，计算首次近似值所使用的常数0x5f3759df也是重要的线索。为确定程序员最初选此常数以近似求取平方根倒数的方法，Charles McEniry首先检验了在代码中选择任意常数R所求取出的首次近似值的精度。回想上一节关于整数和浮点数表示的比较：对于同样的32位二进制数码，若为浮点数表示时实际数值为\textstyle x=(1+m_x)2^{e_x}，而若为整数表示时实际数值则为\textstyle I_x=E_xL+M_x，其中\textstyle 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的相关指数消去，从而将乘方运算化为加法运算。而由于\textstyle \log_2{(x)}与\textstyle \log_2{(x^{-1/2})}线性相关，因此在\textstyle x与\textstyle y_0（即输入值与首次近似值）间就可以线性组合的方式建立方程。在此McEniry再度引入新数\sigma描述\textstyle \log_2{(1+x)}与近似值R间的误差：由于\textstyle 0 \le x ，有\textstyle \log_2{(1+x)}\approx {x}，则在此可定义\sigma与x的关系为\textstyle \log_2{(1+x)}\cong x+\sigma，这一定义就能提供二进制对数的首次精度值（此处0\le\sigma\le\tfrac{1}{3}；当R为0x5f3759df时，有\textstyle \sigma=0.0450461875791687011756）。由此将\textstyle \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.
最后导出的等式\textstyle I_y=R-\frac{1}{2}I_x即与上节代码中一行相契合，由此可见，在平方根倒数速算法中，对浮点数进行一次移位操作与整数减法，就可以可靠地输出一个浮点数的对应近似值。到此为止，McEniry只证明了，在常数R的辅助下，可近似求取浮点数的平方根倒数，但仍未能确定代码中的R值的选取方法。

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

精确度
。]]

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

牛顿法提高精度
在进行了如上的整数操作之后，示例程序再度将被转为长整型的浮点数回转为浮点数（对应），并对其进行一次浮点运算操作（对应），这里的浮点运算操作就是对其进行一次牛顿法迭代，若以此例说明：

: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})，对应的代码即为。

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

历史与考究
的创始人约翰·卡马克。这段代码虽非他所作，但常被认为与他相关。]]

《雷神之锤III》的代码直到QuakeCon 2005才正式放出，但早在2002年（或2003年）时，平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了。

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

注释
参考
脚注
参考文献
*
*
*
*
*
*
*
*

评论 (0)