sqrt函数实现之卡马克方法

sqrt函数的实现主要有三种方式：

1. 二分法
2. 牛顿法
3. 卡马克方法

卡马克方法

这里主要介绍高效的卡马克方法。卡马克方法起源于《雷神之锤III竞技场》中使用的平方根倒数速算法，下列代码是平方根倒数速算法在《雷神之锤III竞技场》源代码中的应用实例。示例剥离了C语言预处理器的指令，但附上了原有的注释：

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与之相减。然而对于代码阅读者来说，他们却难以立即领悟出使用这一常数的目的，因此和其它在代码中出现的难以理解的常数一样，这一常数亦被称为“魔术数字”。如此将浮点数当作整数先位移后减法，所得的浮点数结果即是对输入数字的平方根倒数的粗略估计值，而后再进行一次牛顿迭代法，以使之更精确后，代码即执行完毕。由于算法所生成的用于输入牛顿法的首次近似值已经相当精确，此算法所得近似值的精度已可接受，而若使用与《雷神之锤III竞技场》同为1999年发布的Pentium III中的SSE指令rsqrtss计算，则计算平方根倒数的收敛速度更慢，精度也更低。

将浮点数转化为整数

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

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

"魔术数字"

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

对等式的两边取二进制对数，有

\[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的相关指数消去，从而[将乘方运算化为加法运算。而由于\(log_2(x)\)与\(log_2(x^{-1/2})\)线性相关，因此在\(x\)与\(y_0\)（即输入值与首次近似值）间就可以线性组合的方式创建方程。在此McEniry再度引入新数\(\sigma\)描述\(log_2(1+x)\)与近似值R间的误差：由于\(0≤x<1\)，有\(log_2(1+x)≈x\)，则在此可定义\(\sigma\)与\(x\)的关系为\(log_2(1+x)=x+\sigma\),这一定义就能提供二进制对数的首次精度值（此处\(0≤\sigma≤1/3\)；当R为0x5f3759df时，有\(\sigma=0.0450461875791687011756\)）。由此将\(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，若将其作为长整型表示则示值为\(I_x=E_xL+M_x\)，由此即可根据x的整数表示导出\(y\)（在此\(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. \]

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

牛顿法

在进行了如上的整数操作之后，示例程序再度将被转为长整型的浮点数回转为浮点数（对应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);。

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