快速平方根倒数算法深度理解

快速平方根倒数算法深度理解

快速平方根倒数算法是什么？
简单来说这个算法避开了开方和除法运算快速实现了
$$y=\frac{1}{\sqrt{x}}$$

快速平方根倒数算法首次出现是在著名的雷神之锤3的图形计算优化代码中
话不多说先上c代码
一.代码

小白的我以为这么写的：
float y = 1/sqrt(x);

开发游戏的大佬是这么写的：
float Q_rsqrt(float number)
{
	long i; 
	float x2,y;
	const float threehalfs = 1.5F;
	x2= number * 0.5F;
	y= number;
	i= * (long * ) &y;
	i= 0x5f3759df-(i>>1);
	y= * (float * ) &i;
	//上述代码给牛顿迭代法提供了一个较好的初始值
	y=y * (threehalfs - (x2*y*y)); //牛顿迭代一次   
	//y=y * (threehalfs - (x2*y*y))    可持续迭代
}

经大佬测试代码比(float)(1/sqrt(x))快4倍。

初步观看代码：
①定义了32位整型变量i
②定义了两个32位浮点数
③将1.5存入一个常量
④将形参number除2存入x2
⑤将形参number存入y

之后的四行代码乍一看？？？莫慌下面开始正文

二.理解算法前的知识铺垫
1.IEEE 754 （二进制浮点数算术标准）

32位分为三个部分：符号位 指数部分

①符号位：0为正数 1为负数

②8位指数部分：8位可以表示0到255所有的数字，但我们还需要得到负指数。
因此我们对0-255减去127就可以得到-127到128之间所有的整数。
举个例子如果指数部分我想要4
那这八位必须表示4+127=131 即10000011
这样通过移码我们可以将正负指数都变成0-255之间的整数
③23位尾数：利用23位我们可以表示0到2^23-1之间的整数
对于十进制科学计数法，尾数是1到10
对于二进制科学计数法，尾数是1-2
$$11000=1.1\ast 2^4$$
$$0.0101=1.01\ast 2^{-3}$$
在科学计数法中，第一位按照定义不可以是0，对于二进制这个数只能是1。因此我们不存储这一位。

说了一堆理论的我们举一个例子看一下浮点数是怎么存储在计算机中的：
我们拿 8.25 举例子

Ⅰ.写成二进制：1000.01
Ⅱ.变成科学计数法：$$1.00001\ast 2^3=(1+0.00001)\ast 2^3$$
Ⅲ.解释
指数部分我们用 8bit 表示 3+127
3对应的就是2^3中的3
(1 + 0.00001)：加号前面的1我们丢掉不存
对于0.00001中小数点后的数我们用23位尾数去存

罗嗦了这么久现在我们可以轻松理解以下内容：
对于一个浮点数x(十进制或者二进制)有：
$$x=(-1{)}^s\ast (1+m)\ast 2^e$$

对于下式$$(1+0.00001)\ast 2^3$$
$$m:00001$$
$$e:3$$
我们将存取之后的8位指数部分用E表示
将存取之后的23位尾数部分用M表示

$$E:3+127$$
$$M:00001\cdots （补18个0填满23位）$$
我们可以得出E e和M m的关系式
$$E=e+127$$
$$M=2^{23}\ast m$$
三.数学推导
$$y=\frac{1}{\sqrt{x}}$$
取对数：
$$①{\log }_{2}y=-0.5\ast {\log }_2$$
将y和x分别用二进制科学计数法表示：
$$②y=(1+m_y)\ast 2^{e_y}$$
$$③x=(1+m_x)\ast 2^{e_x}$$

联合上述三个式子化简：
$$lo{g}_2(1+m_y)+e_y=-0.5\ast (lo{g}_2(1+m_x)+e_x)$$
通过一次线性近似：
$$m_y+b+e_y=-0.5\ast (m_x+b+e_x)$$
$$m_y=\frac{M_y}{2^{23}}$$
$$m_x=\frac{M_x}{2^{23}}$$
$$E_y=e_y+127$$
$$E_x=e_x+127$$
上述五个式化简：
$$M_y+2^{23}\ast E_y=1.5\ast 2^{23}\ast (127-b)-0.5\ast (M_x+2^{23}\ast E_x)$$
令
$$I_y=M_y+2^{23}\ast E_y$$
$$I_x=M_x+2^{23}\ast E_x$$
对于$$1.5\ast 2^{23}\ast (127-b)$$我们令其为K，可以知道通过调节b的大小就可以改变K的值

前辈们给出了当
$$b=0.0450465$$时有K=
$$K=0x5f3759df$$
这刚好对应代码中：

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

通过一系列不怎么复杂的数学推导我们究竟得到了什么?
我们知道了一个输出结果的整形和输入整形之间的关系，即：
$$I_y=K-0.5\ast I_x$$
四.代码解释
①

long i; 
float x2,y;
const float threehalfs = 1.5F;
x2= number * 0.5F;
y= number;
i= * (long * ) &y;

我们将数字number存入y中准备对其进行位操作，但是众所周知浮点数没有位操作的语句。
但是长整型数可以进行位操作！

如果执行下面这一行代码：

i=long(y);
y=4.456时i=4

但是我们想要的是将y逐位转化成长整型，换句话说按照上面的代码我们改变了二进制逐位的值。
因此我们采取下面代码：

i= (long * ) &y;

这句话干了什么呢？

简单来说它转换了y对应的内存地址而不是y这个数字本身。
或者说这句话告诉cpu把y指向的内存数据当成长整形来读取

&y              取到浮点数y的地址
(long * ) &y    将该地址转化成长整型
*(long * ) &y    再取其值

6666指针确实是个神奇的东西，我们继续向下
通过这句话浮点数y的二进制表示存入到了i中。
②

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

这句的理解在第三大部分可以轻松找到。
简单来说我们通过放缩系数和平移操作代替了除法和开方的过程，这大大提高了代码运行效率。
③

y = * (float *) &i;

和①类似我们再次通过这种操作将二进制位转换成浮点数。

通过上述三步我们得到了一个近似值。因此我们还需要进行牛顿迭代法继续使结果精确。
④

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

下面进行牛顿迭代法的数学推导：
$$y=\frac{1}{\sqrt{x}}$$
$$f(y)=\frac{1}{y^2}-x$$
$$y_{n+1}=y_n-\frac{\frac{1}{y_n^2}-x}{\frac{-2}{y_n^3}}=y_n\ast (\frac{3}{2}-\frac{1}{2}\ast x\ast y_n^2)$$

五.总结
短短的几行代码背后是强大的数学支撑，万万想不到之前学过的高数会在算法中这样体现。
第一次试写涉及到数学公式推导的文章，如有疑问请大家多多留言！