快速逆平方根
浮点数表示
32位浮点数表示为:
| 符号位 | 阶码 | 尾数 |
|---|---|---|
| s(1) | e(8) | m(23) |
[egin{align}
E=127+e \
M = 10^{23}m \ \
end{align}
]
得到浮点数为：
[x=(-1)^s imes2^{e_x} imes(1+m_x)
]
应用
计算平方根倒数，即：
[y=frac{1}{sqrt{x}}
]
[y^2=frac{1}{x}
]
对最开始的式子两边取对数
[egin{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})也展开，并带入前面浮点数公式得到最终式子：
[egin{align}
m_y+e_y+b=-frac{1}{2}(m_x+b+e_x) \
frac{M_y}{L}+b+E_y-Bapprox{-frac{1}{2}(frac{M_x}{L}+b+E_x-B)} \
M_y+LE_yapprox{frac{3}{2}L(B-b)-frac{1}{2}(M_x+LE_x)}
end{align}
]
所以最终的结果，用(I_y)来表示
[I_yapprox{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;
}