快速平方根倒数算法原理及其动机

原始算法

求某个数的平方根倒数是一个十分常见且高频的操作，比起直接调用数学库里的sqrt函数，一个显著快速的方法是这样的

float InvSqrt(float x){
	float xhalf = 0.5f*x;
	int i = *(int*)&x;
	i = 0x5f3759df - (i>>1);
	x = *(float*)&i;
	x = x*(1.5f-xhalf*x*x); 
	return x;
}

实际上这个算法可以说是为了性能完全牺牲了解释性，因为绝大多数人第一次看到它的时候可能完全摸不着头脑。本文除了详细解释算法的原理以外，还希望揣测算法编写者的动机，并在其他的数值算法性能的的优化上给出指导。

原理

这个算法最关键的两个原理是浮点数的结构以及牛顿法数值迭代的原理。对于前者即算法中的第一步，要理解算法中对浮点数原始二进制数据直接转换为整型后的操作，就是要看对于整型的操作映射到浮点数上的操作是什么，即解释与之间的映射关系。根据IEEE754定义的单精度浮点数的格式应该是这样的

一个单精度浮点数由一位符号位，八位指数位以及二十三位小数位组成，若将一位符号位，指数位作为有符号8位整型，小数位作为无符号23位整型，并分别以表示，则对应表示的单精度浮点数为

上述的符号位在这个算法内总是代表正数，因此可以忽略。如果将上述原始二进制直接转换为32位整型（忽略了浮点数的符号位转换后实际上是一个31位无符号数），则其表示为

观察到中小数位和指数位是相乘的形式，而中则是乘以一个系数后相加的线性形式，很容易联想到它们之间的映射应该是指数/对数的映射。对两边求2的对数可得

同时考虑到较为接近0，可以用泰勒展开的一阶近似。因此根据上式可以大致认为浮点数的表示和它原始数据所直接对应的整型是一个对数映射关系。
采用对数映射后得到的整型恰好有一个非常优越的性质，在计算平方根时可以直接对整型除以2得到其大致的在整型上的结果。但是目前这种映射实际上采用了过多近似，其精度不甚理想。在原有方法的基础上，我们可以把在原点的泰勒展开的近似增加一个偏置，在整个0与1之间的范围上精度能够更好。
将上述的一阶近似改写为，其中是一个可优化的参数，在不同判定拟合精度的指标（例如平方和最小，最大误差等）下的最优值，则此时

接下来如果令平方根倒数的估计值为则，由上述转换为整型的近似关系可知

当时，第一项，而后面一项除以2则可以使用移位操作完成，将原始数据直接转换回浮点数即可得到一个估计的平方根倒数。

在已知上述估计值的情况下，进一步提高精度就要进行第二步，这里就是数值迭代方法中的经典方法牛顿法了。令，则的根即为，牛顿法的迭代式可以表示为

对应代码内返回前一行的迭代计算。实际上由于上述估计值已经比较接近真实值，牛顿法在此处收敛较快，在这里迭代一次的精度就已经足够了。

总的来说，计算平方根倒数首先是利用直接转换整型，并进行计算的操作的得到一个平方根倒数的估计值，然后再利用牛顿法迭代一次确保最终的精度。这个算法的显著特征就是速度极快，另外还需要优化的就是其精度，例如的选值，牛顿法迭代的求根函数等。还有一种玄学调参的思维是，将上述算法中的常量全部作为优化的变量，暴力搜索寻找一组使得误差在一个范围内总体最小的参数，这里有人做过。但是据说原始的0x5f3759df这个数目前没有办法解释如何找到它，按照其他优化指标计算出的常量算法的精度都比原始值优秀。

动机

上面提到的是对算法本身原理的逆向讨论，即已知有人写出过这样的算法我们如何解释这个算法的原理。但更重要的了解编写者的动机，即怎么想到一个这样的算法能够极大地优化平方根倒数的计算。
在我看来，一个优化的关键是优化是借助于计算机浮点数结构完成的。同一段32位数据用浮点数表示和用整型表示，在我们看来数的意义是不一致的，而对于计算机来说是完全没有区别的。因此由于浮点数指数项的存在，借助这种转换关系可以完成一种对数/指数的映射。第二个优化的关键思路是在牛顿迭代法前，先用某种更加高效的手段估计一个不那么精确的值作为初始解，或者说牛顿法可以与其他方法组合使用。牛顿迭代法在根附近呈平方收敛，一个合适的初始值可以代替多次迭代达到相同的精度。

快速平方根倒数算法的类似计算方法同样可以计算平方根以及2n次方根等，另外整型的直接转换也为也是一种优化速度的新思路，如果针对定点计算以及要求通用性则可以采用CORDIC计算

参考

Fast inverse square root - wikipedia
IEEE 754 - wikipedia
CHRIS LOMONT - FAST INVERSE SQUARE ROOT