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