快速逆平方根
快速逆平方根
浮点数表示
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;
}
