[算法]快速平方根倒数算法学习-Fast Inverse Square Root

简介：平方根倒数速算法是一种用于快速计算32位浮点数的算术平方根的倒数的算法。该算法最早可能由SGI在90年代前期发明，并在1999年的《雷神之锤III竞技场》源代码中应用，但直到2002-2003年才在Usenet等公共论坛上出现。该算法的优势在于减少了求平方根倒数时浮点运算操作带来的大量运算耗费。在计算机图形学领域，例如计算照明和投影的波动角度与反射效果时，常需要计算平方根倒数。算法通过对32位带符浮点数进行一系列操作，包括右移、使用特定十六进制魔术数字的减法以及牛顿迭代法，以快速得到平方根倒数的近似值。该算法在相同精度下比直接使用浮点数除法快四倍。       

先给出平方根倒数算法的代码：

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));   // 2st iteration
										//this can be removed

	return y;
}

到底是什么样的代码能让3D游戏之父、FPS游戏之父、程序员之神约翰·卡马克（John Carmack）直呼“what the fuck?”

在分享这个算法之前，先来算一下根号2的近似值是多少，相比各位都可以秒答是1.414，但是，如果要求近似值精确到后10位，后15位小数呢？这样就有点困难了吧。

不妨先换个思路，求2的平方根，其实就是求y=x2-2(x平方-2)于x轴的交点，这么求呢，先从(2,0)作垂线，找到和y=x2-2的交点，作切线，然后求导得到这条切线的斜率，这样就能很快求出该切线与x轴的交点。再向上作垂线，重复几次，就会发现这个近似值会不断变得越来越精确直到满足要求，而这就是300年前牛顿的方法，把初始值带入，通过作切线，不断迭代，而且因为每一步的方法都一样，所以只有一个公式

用人工计算或许比较麻烦，但是特别适合用计算机来计算：

//牛顿法求平方根
//传入原数和要精确的位数，比如精确到小数后5位就输入1e-5
float newton_sqrt(float number,float precision_) {
	float x = number;
	float precision = precision_;//要求的精度

	while (true) {
		float next_x = 0.5*(x + 2 / x);//牛顿法
		if (abs(next_x - x) < precision) {//判断精度是否达到要求
			return next_x;
		}
		x = next_x;
	}

}

 这种方法只要会写for循环就能实现，但是还有天才想着要简化算法，不用循环，只要一步牛顿法就能得出结果，说回牛顿法，如果求2的根号近似值时，第一步带入近似值1.414而不是2的话，那么只要一步牛顿法即可得到一个非常精确的近似值，也就是说，简化牛顿法的关键在于找到一个非常接近近似值的数字来作为初始值，但是问题就来了，1.414是我本来就知道的，但是在绝大部分情况下，我是因为不知道解才需要用牛顿法来求解，但是你却要我先给出一个这非常接近近似值的数字，不是形成了一个死循环了么！？这就引出了本文的主题，平方根倒数快速算法。

 一个32位的浮点数，第1位是符号位（S），0代表正数，1代表负数。第2-9位是指数位（E），相当于科学计数法的指数，比如5×106的6。第10-32位是有效数字部分（M），相当于科学计数法前面的有效部分。我这里讲的比较简略，如果有不明白计算机是这么存储浮点数的请先去其他地方了解一下先。

总之，一个浮点数的值可以表示为

这么麻烦的存储方式，还牛顿法求倒数，那还不然老老实实用传统牛顿法，但是程序员大佬Greg.Walsh格雷格·沃什在解决这个问题的时候引入了对数，，而对数的运算强在能把乘除运算转化位加减运算，

  

到这里，我们基本能探明Greg.Walsh这四行代码的想法：

float Q_rsqrt(float number) {//快速平方根倒数算法
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y = number;

	i = * (long * ) &y; //把y强转为可以进行二进制运算的long类型
	i = 0x5f3759df - (i >> 1); //计算y平方根倒数的近似值i
	y = * (float * ) &i; //把i强转为浮点数
	y = y * (threehalfs - (x2 * y * y)); //使用牛顿法计算精确值
	//y = y * (threehalfs - (x2 * y * y)); //一步就已经足够精确，不用第二步
	return y;
}

而这个算法之后经常被用于3D游戏当中，因为再计算机各种模型曲面，本质上都是海量的三角面，要渲染光照再物体表面反射的的效果，那就要知道这个三角面的法向量，法向量的长度平方R很容易被求出，则单位法向量=坐标值/R的平方根。

 其他话：这是我再CSDN的第四篇文章，也是第一篇讲解算法的文章，CSDN里面自带的公式模板等还没有使用得太明白，所以很多文本都是用了截图，请见谅。这个快速平方根倒数算法也是我在看了其他视频、文章后，特别是什么代码让程序员之神感叹“卧槽”？改变游戏行业的平方根倒数算法_哔哩哔哩_bilibili

经过自己的学习理解后再用自己的话语讲解出来，不知道有没有一些拾人牙慧的嫌疑，不过教会他人是学习的最好方法，为了教会，必须对已有的知识进行层层梳理、理解，这次分享也是我对于这个算法有了深刻理解，讲得不好的地方还请读者在评论区讨论交流。