快速逆平方根
浮点数表示
32位浮点数表示为:
符号位 | 阶码 | 尾数 |
---|---|---|
s(1) | e(8) | m(23) |
得到浮点数为:
\[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;
}