快速逆平方根

浮点数表示

32位浮点数表示为:

符号位	阶码	尾数
s(1)	e(8)	m(23)

\[\begin{align} E=127+e \\ M = 10^{23}m \\ \\ \end{align} \]

得到浮点数为：

\[x=(-1)^s\times2^{e_x}\times(1+m_x) \]

应用

计算平方根倒数，即：

\[y=\frac{1}{\sqrt{x}} \]

\[y^2=\frac{1}{x} \]

对最开始的式子两边取对数

\[\begin{align} \log_{2}{y}=-\frac{1}{2}\log_{2}{x} \\ \log_{2}{y}=-\frac{1}{2}[\log_{2}{(1+m_x)}+e_x] \end{align} \]

对于\(y=\log_{2}{(1+m_x)}\)​​可以画出一个平滑的函数图像：

因为 \(m_x\)​取值范围在\((0,1)\)内

所以可以认为

\[y=\log_{2}{(1+m_x)}\approx{mx+b} \]

这时需要计算出b的值来使得方差最小

计算出来的b的值为：0.0430357

所以原方程变为：

\[-\frac{1}{2}\log_{2}{(1+m_x)}+e=-\frac{1}{2}(m_x+b+e_x) \]

上述中左边\(\log_{2}{y}\)也展开，并带入前面浮点数公式得到最终式子：

\[\begin{align} m_y+e_y+b=-\frac{1}{2}(m_x+b+e_x) \\ \frac{M_y}{L}+b+E_y-B\approx{-\frac{1}{2}(\frac{M_x}{L}+b+E_x-B)} \\ M_y+LE_y\approx{\frac{3}{2}L(B-b)-\frac{1}{2}(M_x+LE_x)} \end{align} \]

所以最终的结果，用\(I_y\)来表示

\[I_y\approx{\frac{3}{2}L(B-b)-\frac{1}{2}I_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;  // 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 ) ); // 2nd iteration, this can be removed
  	return y;
}

核心部分即为代码中第8行表示，也是著名的注释的由来

精简版：

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

Reference

计算机与数学 —— 雷神之锤3源码中的快速逆平方根算法

20年前的几行代码竟如此牛逼？惊了

标签：approx,frac,log,align,float,平方根,快速,浮点数
来源： https://www.cnblogs.com/JoshuaYu/p/15091492.html