6
求平方根倒数的算法
source link: https://zhiqiang.org/coding/inverse-square-root-algorithm-analysis.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
求平方根倒数的算法
下面这个求1/√x1/x 的函数号称比直接调用 sqrt 库函数快 4 倍，来自游戏 Quake III 的源代码。
float InvSqrt (float x){
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i>>1);
float y = *(float*)&i;
y = y*(1.5f - xhalf*y*y);
return y;
}
我们这里分析一下它的原理（指程序的正确性，而不是解释为何快）。
分析程序之前，我们必须解释一下 float 数据在计算机里的表示方式。一般而言，一个 float 数据xx 共 32 个 bit ，和 int 数据一样。其中前 23 位为有效数字MxMx ，后面接着一个 8 位数据ExEx 表示指数，最后一位表示符号，由于这里被开方的数总是大于 0 ，所以我们暂不考虑最后一个符号位。此时
x=1.Mx2Ex−127(1)(1)x=1.Mx2Ex−127
如果我们把计算机内的浮点数xx 看做一个整数IxIx ，那么
Ix=223Ex+Mx(2)(2)Ix=223Ex+Mx
现在开始逐步分析函数。这个函数的主体有四个语句，分别的功能是：
int i = *(int*)&x;这条语句把xx 转成i=Ixi=Ix 。i = 0x5f3759df - (i>>1);这条语句从IxIx 计算I1/√xI1/x 。y = *(float*)&i;这条语句将I1/√xI1/x 转换为1/√x1/x 。y = y\*(1.5f - xhalf\*y\*y);这时候的 y 是近似解；此步就是经典的牛顿迭代法。迭代次数越多越准确。
关键是第二步 i = 0x5f3759df - (i>>1); 这条语句从IxIx 计算I1/√xI1/x ，原理:
令y=1/√xy=1/x ，用x=(1+mx)2exx=(1+mx)2ex 和y=(1+my)2eyy=(1+my)2ey 带入之后两边取对数，再利用近似表示log2(1+z)∼z+δlog2(1+z)∼z+δ ，算一算就得到
Iy=32(127−δ)223−Ix/2(3)(3)Iy=32(127−δ)223−Ix/2
若取δ=0.0450465679168701171875δ=0.0450465679168701171875 ，32(127−δ)22332(127−δ)223 就是程序里所用的常量 0x5f3759df。至于为何选择这个δδ ，则应该是曲线拟合实验的结果。
Q. E. D.
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK