快速逆平方根

快速逆平方根

浮点数表示

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年前的几行代码竟如此牛逼?惊了

上一篇:【剑指offer】65. 不用加减乘除做加法


下一篇:筛选数组元素方法