平方根倒数速算法(Fast inverse square root)，经常和一个十六进制的常量 0x5f3759df联系起来。该算法被用来快速运算平方根倒数，速度是 float(1/sqrt(x)) 方法的4倍。该算法大概由上个世纪90年代的硅图公司开发出来，后来出现在John Carmark的Quake III Arena的源码中。
这是一个古老的算法，最早的讨论见于2001年中国的CSDN论坛上。并且该段代码可能已经不适用于当代的64bits机器，因为现在的64bits的机器上 long 型数据长度已经从4字节变成了8字节。但是我觉得，一代代新人对于一些老的经典的问题的讨论是有实际应用价值以外的意义的。仅此而已，不做争辩。
代码总览
下面的快速平方根倒数代码来自Quake III Arena。
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展示了其在其他浮点存储标准下的运用。
算法运行示例
这里有一个例子，令x = 0.15625，我们想计算的值
。程序的第一步如下：
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??
算法原理
该算法通过下面几个步骤计算
的值：
- 将32bits的浮点型变量reinterpret成整型变量，用于计算log2(x)的近似值
- 用上述的近似值计算log2(1/√x)的近似值
- 将上述的计算结果强转回32bits浮点型变量，得到迭代初值
- 经过一次简单的牛顿迭代，得到足够高精度的结果
- 将浮点数型数据强制转换成整型，近似求解对数值
前面一篇博文专门介绍了IEEE 754标准下的浮点数表示方法，下面简单的回顾一下。为了能够存储非零的浮点数，第一步就是将浮点数x写成标准化的二进制形式：
这里ex是整数，mx ∈ [0, 1)显而易见。我们把1.b1b2b3...称为尾数1 + mx的二进制表示，其中小数点左边的1称为前导数位。很显然，前导数位始终为1，所以在浮点数的存储过程中将其忽略。
我们现在来计算三个数：
- Sx，符号位。当x>0，等于0，x<0，等于1；
- Ex = ex + B，偏移指数。为了既能存储正指数，又能存储负指数。IEEE 754标准在真实的指数上加一个偏移量来作为存储指数，对于单精度浮点数来说B=127；
- Mx = mx × L，L = 223
首先一个被开方数肯定是个非负数，所以符号位一定为0，那么
，对其取以2为底的对数可得
。因为mx ∈ [0, 1)，所以有
。其中σ是用来调节近似的精确度的，当σ ≈ 0.0430357时可以获得最优的近似结果。所以我们有
如果我们把浮点数强制转换成整数的话，偏移指数变成了这个整数的高8位，尾数部分变成了整数的低23位。
这个整数的数值为：
从上式我们可以看出，整数值Ix其实是对log2(x)平移和伸缩后的结果，并且很容易把log2(x)重新表示出来：
- 迭代初值的计算
好的，到现在为止我们知道了将浮点数强转成整数发生了什么。卡马克开方法最重要的一步就是找出了一个迭代初值，这个初值神接近真实结果，使得迭代收敛速度奇快。那么我们接下来将要看到，所谓的Magic Number的来历，以及这个迭代初值是怎么求得的。
令y = 1/√x ，然后就有
根据上面刚刚得到的结论有：
整理后得到：
我们发现了什么？将平方根倒数和被开方数强转成整数竟然有这种近似关系！所以所谓的Magic Number就是3/2L(B-σ)的十六进制写法，所以就有了如下的代码：
i = 0x5f3759df - ( i >> 1 );
第二项是通过右移操作计算½ Ix。
- 牛顿迭代法
通过上述的操作后，该算法将求得的整数重新强转回浮点数，这样就得到了迭代的初始值。下面再来看看迭代是如何进行的。
令
y就是就是下面方程的解：
牛顿迭代法给出了下面的方法来迭代求解方程的根：
其中
是前一次迭代的结果，对于第一次迭代来说就是迭代初值。所以上式写成代码的形式如下：
y = y*(threehalfs - x2*y*y);
重复上面的步骤，将每次计算的结果
作为输入，就会不断逼近真实结果。
由于我们通过前面的方法选取的初始值已经特别接近真实值了，没有必要进行更多次的迭代，因此源码中甚至将第二次迭代的过程都注释掉了。
Reference:
[1] https://en.wikipedia.org/wiki/Fast_inverse_square_root