平方根倒数速算法

要理解这段代码，首先需了解浮点数的存储格式。一个浮点数以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，而后面的各位则用于表示其数值。将浮点数取别名存储为整数时，该整数的数值即为${\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]。

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

如上所述，平方根倒数速算法所得的近似值惊人的精确，右图亦展示了以上述代码计算（以平方根倒数速算法计算后再进行一次牛顿法迭代）所得近似值的误差：当输入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]。