sqrt函数实现之卡马克方法

Posted 2021-01-06 muyangshaonian

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位表示的则是有效数字中除最高位以外的其余数字。将上述结构表示成公式即为，其中{displaystyle scriptstyle m}表示有效数字的尾数（此处{displaystyle scriptstyle 0leq m<1}，偏移量{displaystyle scriptstyle B=127}[文 10]，而指数的值{displaystyle scriptstyle E-B}决定了有效数字（在Lomont和McEniry的论文中称为“尾数”（mantissa））代表的是小数还是整数[文 11]。以上图为例，将描述带入有{displaystyle scriptstyle m=1 imes 2^{-2}=0.250}），且{displaystyle scriptstyle E-B=124-127=-3}，则可得其表示的浮点数为{displaystyle scriptstyle x=(1+0.250)cdot 2^{-3}=0.15625}。