平方根倒数算法
平方根倒数算法
看到一个非常巧妙的平方根倒数计算算法，特此记录。
我们首先要先了解IEEE754浮点数标准：32位的浮点数，划分为符号位，阶码，尾数。我们将一个小数以二进制的科学记数法表示，如：1.1 * 2^3。阶码即指数部分，尾数即小数点后面的数字（小数点前面一定是1：二进制小数点前面只能是0/1，同时科学计数法又要求该数大于等于1，因此可以省略）。
- 符号位占1位，0/1表示该数的正/负
- 阶码占8位，由于需要表示负指数，因此范围为[-127, 128]
- 尾数占23位
我们将阶码记作E，尾数记作M，由此得到公式：
$$
(1+\frac{M}{2^{23}})*2^{E-127}
$$
接下来说回求平方根算法。早期的平方根算法利用了牛顿求解零点的迭代法。其原理是先任取一点x0，然后计算x0该点的切线与x轴的交点，将该交点作为下一次迭代的x0。由此迭代，x0将逼近函数与x轴的真实交点。形式上，迭代公式为：
$$
x_{n+1}=x_{n}-\frac{f(x_{n})}{f’(x_n)}
$$
最原始的计算算法需要多次迭代，直至精度满足需求。但是一些聪明的程序员发现初始值的选取很重要：如果我们要计算根号2，初始值选了2，需要迭代五六次才能达到1e-7的精度需求；而如果我们初始值选了1.414，一次迭代便可达到1e-7的精度需求。那么现在的问题就变成了：如果找到一个比较近似的解？
我们要计算$$a^{-\frac{1}{2}}$$，考虑它的对数，即为$$-\frac{1}{2}*log_2{a}$$。那么我们如何求$$log_2{a}$$呢，代入上面计算机存储浮点数的公式，我们需要求的即：
$$
log_2({(1+\frac{M}{2^{23}})*2^{E-127}})=log_2{(1+\frac{M}{2^{23}})}+E-127
$$
那么如何求$$log_2{(1+\frac{M}{2^{23}})}$$呢？考虑到$$\frac{M}{2^{23}}$$是个介于（0，1）的数，而我们注意到y=log_2(1+x)的函数图像在（0，1）区间与y=x非常接近，这样我们不如直接取$$\frac{M}{2^{23}}$$作为其近似值。由此我们得到初始值的选取公式：
$$
\frac{M}{2^{23}}+E-127=\frac{1}{2^{23}}(2^{23}E+M)-127
$$
注意到，E2^23正好是阶码左移23位置，再加上M即浮点数在内存中的表示。我们设其内存中表示的这个数字为Y，原浮点数为a，那么就有：
$$
log_2{a}=\frac{1}{2^{23}}*Y-127
$$
代入上面，有：
$$
log_2{a^{-\frac{1}{2}}}=-\frac{1}{2}(\frac{1}{2^{23}}*Y-127)
$$
我们设$$a^{-\frac{1}{2}}$$这个浮点数在内存中表示为A，那么有：
$$
\frac{1}{2^{23}}A-127=-\frac{1}{2}(\frac{1}{2^{23}}Y-127)
$$
解得：
$$
A=3812^{22}-\frac{1}{2}Y
$$
3812^22在16进制的表示下为0x5f400000，注意到我们之前用y=x近似对数，实际上可以将y=x向上平移一些，取得更加精确的结果，最终这个值定为0x5f3759df。
那么最终的程序如下：
1 | float Q_rsqrt(float number){ |
