平方根算法
希望你认真看完这篇文章，我想一定会启发到你
引言：QuakeIII 的平方根倒数
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;
}
这恐怕是程序员界里最出名的算法之一，它最奇特的地方在于使用了一个 Magic Number 即 0x5f3759df 但是，这个数字是怎么来的呢，这个算法有哪些可以应用到我们的平方根算法中的思路呢？而这就是我们要讨论的。
本文章需要你至少理解“求导”和“牛顿迭代”，我不会在这里赘述，如果不会，自行百度。
计算机中的浮点数
理解上面算法的前提是明白计算机如何存储浮点数，请自学 ieee754 浮点数标准。
注意：尾数位的数字是大于等于0小于等于1的
这两个公式就是算法操作的基础。
指数偏移量 = 2 指数位位数 − 1 指数偏移量 = 2^{指数位位数} - 1 指数偏移量=2指数位位数−1
浮点数 = − 1 符号位 × 2 指数位 − 指数偏移量 × ( 1 + ∑ i = 0 m = 尾数位比特数 m m − i 2 − i ) 浮点数 = -1^{符号位}\times2^{指数位-指数偏移量}\times (1+\sum_{i=0}^{m=尾数位比特数}{m_{m-i}2^{-i}} ) 浮点数=−1符号位×2指数位−指数偏移量×(1+∑i=0m=尾数位比特数mm−i2−i)
这里我们仅针对32位浮点数（64位是类似的）并且令指数位为 E E E 尾数位为 M M M并且由于是平方根，所以符号位总是0，以上的公式就可以化为下面的:
浮点数 = 2 E − 127 × ( 1 + M ) 浮点数 = 2^{E-127}\times(1+M) 浮点数=2E−127×(1+M)
算法分析
估数
y = 1 x = x − 1 2 x = 2 E − 127 ∗ ( 1 + M ) 所以 y = 2 − 1 2 ( E − 127 ) ∗ ( 1 + M ) − 1 2 令 E ′ = y 的指数位 令 M ′ = y 的尾数位 那么 y = 2 E ′ − 127 ∗ ( 1 + M ′ ) 因为 0 < = M < = 1 并且 0 < = M ′ < = 1 所以 1 2 < = ( 1 + M ) − 1 2 < = 1 并且 1 < = ( 1 + M ′ ) < = 2 若是不考虑尾数，只考虑指数我们就可以近似的得到 2 E ′ − 127 ≈ 2 − 1 2 ( E − 127 ) 即 E ′ ≈ 190.5 − 1 2 E y = \frac{1}{\sqrt{x}}=x^{-\frac{1}{2}}\\ x = 2^{E-127}*(1+M)\\ 所以 y = 2^{-\frac{1}{2}(E-127)}*(1+M)^{-\frac{1}{2}}\\ 令E' = y的指数位\\ 令M' = y的尾数位\\ 那么 y = 2^{E'-127}*(1+M')\\ 因为 0<=M<=1并且0<=M'<=1\\ 所以 \frac{1}{\sqrt{2}}<=(1+M)^{-\frac{1}{2}}<=1\\ 并且 1<=(1+M')<=2\\ 若是不考虑尾数，只考虑指数我们就可以近似的得到\\ 2^{E'-127}\approx2^{-\frac{1}{2}(E-127)}\\ 即E'\approx190.5-\frac{1}{2}E y=x1=x−21x=2E−127∗(1+M)所以y=2−21(E−127)∗(1+M)−21令E′=y的指数位令M′=y的尾数位那么y=2E′−127∗(1+M′)因为0<=M<=1并且0<=M′<=1所以21<=(1+M)−21<=1并且1<=(1+M′)<=2若是不考虑尾数，只考虑指数我们就可以近似的得到2E′−127≈2−21(E−127)即E′≈190.5−21E
$$
得到了
E
E
E和
E
′
E'
E′的关系后，我们就可以开始神奇的魔法，但是由于我们的数字是浮点数，所以需要先将其转化为整数，以便后续可以进行位的操作，而这就是代码中i = * ( long * ) &y;做的是，就是用整数形式来表示这个浮点数。
注意我所说的形式因为往后的所有操作都是形式上的
如此一来我们便可以将190表示为二进制数10111110将其作为一个常数放在指数为，就能得到y' = 0x5f000000-(x>>1);就是
E
′
E'
E′的整数形式，这个形式也是32位的并且在 hack时候记得要在数字前补一个0，作为符号位。
可以看到0x5f000000已经很接近那个 Magic Number了。
0x5f3759df这个数字剩余的部分，实际就是对尾数位进行了一些偏移，以进一步拉近精度，这里就不再进行详细推导了，可以通过改变值，多试几次，试出一个还不错的。
精确
以上的步骤得到的 y 其实已经和实际的结果有将近 60%-70% 的接近了。那么最后就是再将整数形式的浮点数换回真正的浮点数，然后使用牛顿迭代进一步逼近值。
对于
f
(
x
)
=
1
x
2
−
n
f(x) = \frac{1}{x^2}-n
f(x)=x21−n它的根实际就是
1
n
\frac{1}{\sqrt{n}}
n1而这个函数的导函数就是
f
′
(
x
)
=
−
2
x
−
3
f'(x)=-2x^{-3}
f′(x)=−2x−3
根据牛顿迭代公式 x n + 1 = x n − f ( x ) f ′ ( x ) x_{n+1}=x_n-\frac{f(x)}{f'(x)} xn+1=xn−f′(x)f(x)这个式子可以变形为 x n + 1 = 1.5 ( x n + n ∗ x n 3 ) x_{n+1}=1.5(x_n+n*x_n^3) xn+1=1.5(xn+n∗xn3)
于是我们就可以进一步精确数字，上面那个算法的思路也就完整剖析出来了。
实现
仿照以上的思路，我们也可以完成我们的平方根函数，思路是类似的。如果你还能看到这里，那么我鼓励你自己拿起笔计算一下，我不会再计算了，直接给出代码。对于编程来说，始终热爱思考是十分重要的。
// 仿 Quake 的浮点 hack 平方根。
float Q_sqrt(float n) {
float y = n;
long i = *(long*)&y;
i = 0x1fbb4eee + (i >> 1); // 对浮点整形表示进行hack
y = *(float*)&i; // 再转回浮点型
y = 0.5 * (y + n / y); // 一次迭代
y = 0.5 * (y + n / y); // 找的数误差还有点大，所以进行第二次迭代...
// 对 0 进行额外检查，不然它会得到一个 1.xxe-20 这样一个非常接近0的小数。
return n == 0 ? 0 : y;
}
这样，算法保证速度，而且精度能保证小数点后4位都是准确的，但是精度和速度都仍然可以在此之上继续优化。
参考 ： 没那么神秘的快速平方根倒数算法
扫一扫
2373
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
