sqrt 函数的实现主要有三种方式:
二分法
牛顿法
卡马克方法
这里主要介绍高效的卡马克方法. 卡马克方法起源于《雷神之锤 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 0\leq m<1}, 偏移量{\displaystyle \scriptstyle B=127}[文 10], 而指数的值{\displaystyle \scriptstyle E-B} 决定了有效数字 (在 Lomont 和 McEniry 的论文中称为 "尾数"(mantissa)) 代表的是小数还是整数[文 11]. 以上图为例, 将描述带入有{\displaystyle \scriptstyle m=1\times 2^{-2}=0.250}), 且{\displaystyle \scriptstyle E-B=124-127=-3}, 则可得其表示的浮点数为{\displaystyle \scriptstyle x=(1+0.250)\cdot 2^{-3}=0.15625}.
来源: http://www.bubuko.com/infodetail-2768428.html