平方根倒数速算法（卡马克开方法）

时间 2019-11-07

标签平方根倒数速算马克开方繁體版

原文原文链接

平方根倒数速算法(Fast inverse square root)，常常和一个十六进制的常量 0x5f3759df联系起来。该算法被用来快速运算平方根倒数，速度是 float(1/sqrt(x)) 方法的4倍。该算法大概由上个世纪90年代的硅图公司开发出来，后来出如今John Carmark的Quake III Arena的源码中。html

这是一个古老的算法，最先的讨论见于2001年中国的CSDN论坛上。而且该段代码可能已经不适用于当代的64bits机器，由于如今的64bits的机器上 long 型数据长度已经从4字节变成了8字节。可是我以为，一代代新人对于一些老的经典的问题的讨论是有实际应用价值之外的意义的。仅此而已，不作争辩。算法

代码总览

下面的快速平方根倒数代码来自Quake III Arena。this

 1 float Q_rsqrt( float number )
 2 {
 3     long i;
 4     float x2, y;
 5     const float threehalfs = 1.5F;
 6 
 7     x2 = number * 0.5F;
 8     y  = number;
 9     i  = * ( long * ) &y;                       // evil floating point bit level hacking
10     i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
11     y  = * ( float * ) &i;
12     y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
13 //    y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
14 
15     return y;
16 }

为了求解平方根的倒数值，首先须要产生一个近似值。而后在利用数值计算的方法来减少近似值的偏差。近似值的选取很大程度上影响后面数值计算迭代次数。在这以前，这个近似值的选取是经过查表获得的，利用所谓的lookup table。但时上面提供了比查表还要迅速的方法，而且速度是常规计算平方根倒数的四倍。虽然该方法在精度上有必定的损失，但时计算速度上弥补了这一点的不足。该算法是为IEEE 754标准下的浮点数存储规则设计的，Chris Lomot和Charles McEniry展现了其在其余浮点存储标准下的运用。spa

算法运行示例

这里有一个例子，令x = 0.15625，咱们想计算的值。程序的第一步以下：.net

0011_1110_0010_0000_0000_0000_0000_0000  //x和i的各比特位表示
0001_1111_0001_0000_0000_0000_0000_0000  //右移一位后的结果: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  //magic number 0x5f3759df
0100_0000_0010_0111_0101_1001_1101_1111  //0x5f3759df - (i >> 1)的结果

将上述的数若是是用IEEE 754标准存储的浮点数的话，各个数值大小以下：设计

0_01111100_01000000000000000000000  //1.25 * 2^-3
0_00111110_00100000000000000000000  //1.125 * 2^-65
0_10111110_01101110101100111011111  //1.432430... * 2^+63
0_10000000_01001110101100111011111  //1.307430... * 2^+1

咱们神奇的发现最后的数字已经很接近准确值2.52982了！what the fuck, why??3d

算法原理

该算法经过下面几个步骤计算的值：code

将32bits的浮点型变量reinterpret成整型变量，用于计算log₂(x)的近似值
用上述的近似值计算log₂(1/√x)的近似值
将上述的计算结果强转回32bits浮点型变量，获得迭代初值
通过一次简单的牛顿迭代，获得足够高精度的结果

将浮点数型数据强制转换成整型，近似求解对数值

前面一篇博文专门介绍了IEEE 754标准下的浮点数表示方法，下面简单的回顾一下。为了可以存储非零的浮点数，第一步就是将浮点数x写成标准化的二进制形式：htm

${\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}$

这里 $e x 是整数， m x \in [0, 1) 显而易见。咱们把 1.b 1 b 2 b 3 ... 称为尾数 1 + m x 的二进制表示，其中小数点左边的1称为前导数位。很显然，前导数位始终为1，因此在浮点数的存储过程当中将其忽略。$ blog

咱们如今来计算三个数：

$S x ，符号位。当 x>0 ，等于0， x<0 ，等于1；$
$E x = e x + B ，偏移指数。为了既能存储正指数，又能存储负指数。IEEE 754标准在真实的指数上加一个偏移量来做为存储指数，对于单精度浮点数来讲B=127；$
$M x = m x \times L ， L = 2 23$

$首先一个被开方数确定是个非负数，因此符号位必定为0，那么，对其取以2为底的对数可得。由于 m x \in [0, 1)，因此有。其中 σ 是用来调节近似的精确度的，当 σ \approx 0.0430357 时能够得到最优的近似结果。因此咱们有$

$若是咱们把浮点数强制转换成整数的话，偏移指数变成了这个整数的高8位，尾数部分变成了整数的低23位。$

$这个整数的数值为：$

从上式咱们能够看出，整数值 $I x 实际上是对 log2(x) 平移和伸缩后的结果，而且很容易把 log2(x) 从新表示出来：$

迭代初值的计算

好的，到如今为止咱们知道了将浮点数强转成整数发生了什么。卡马克开方法最重要的一步就是找出了一个迭代初值，这个初值神接近真实结果，使得迭代收敛速度奇快。那么咱们接下来将要看到，所谓的Magic Number的来历，以及这个迭代初值是怎么求得的。

令 $y = 1/ \sqrt x$ ，而后就有

根据上面刚刚获得的结论有：

整理后获得：

咱们发现了什么？将平方根倒数和被开方数强转成整数居然有这种近似关系！因此所谓的Magic Number就是3/2L(B-σ)的十六进制写法，因此就有了以下的代码：

i  = 0x5f3759df - ( i >> 1 );

第二项是经过右移操做计算½ I_x_。

牛顿迭代法

经过上述的操做后，该算法将求得的整数从新强转回浮点数，这样就获得了迭代的初始值。下面再来看看迭代是如何进行的。

令

y就是下面方程的解：

牛顿迭代法给出了下面的方法来迭代求解方程的根：

其中 $\,y_{n}$ 是前一次迭代的结果，对于第一次迭代来讲 $\,y_{n}$ 就是迭代初值。因此上式写成代码的形式以下：

y = y*(threehalfs - x2*y*y);

重复上面的步骤，将每次计算的结果做为 $\,y_{n}$ 输入，就会不断逼近真实结果。

因为咱们经过前面的方法选取的初始值已经特别接近真实值了，没有必要进行更屡次的迭代，所以源码中甚至将第二次迭代的过程都注释掉了。

Reference:

[1] https://en.wikipedia.org/wiki/Fast_inverse_square_root