详细解析Fast Inverse Square Root

详细解析Fast Inverse Square Root算法（

目前已被math库收录
学习价值更重要
该内容仅针对32位编码架构

Fast Inverse Square Root算法是用来快速计算如开方的倒数计算：

$$FISqrt(x)=\frac{1}{\sqrt{x}}$$

$$⟺FISqrt(x)=x^{-\frac{1}{2}}$$

零、Fast Inverse Square Root算法代码：

float FISqrt(float const number)
{
	long i;
	float x2;
	float y;
	const float threehalfs = 1.5F;
	x2 = number * 0.5F;
	y = number;
	i = * ( long * ) &y;
	i = 0x5f3759df - ( i >> 1 );
	y = * ( float * ) &i;
	y = y * ( threehalfs - ( x2 * y * y ) );
	return y
}

---

一、背景

（一） IEEE 754 浮点数表示法（数值转二进制码）

要看懂这个代码，就必须了解浮点数是如何在2进制中存储的。IEEE754规定了现代计算机浮点数的存储方式。

IEEE 754 浮点数表示法其实就是二进制版本的科学计数法。

例：

$$1678.4125$$

1.使用十进制科学计数法表示：

$$1.6784125\times 10^3$$

2.使用二进制科学计数法表示：

（1）将 $1678.4125$ 转换为二进制数：

    整数部分：$1678_{(10)}=11010001110_{(2)}$

    为了学习我们进行手动计算：

    小数部分：$0.4125_{(10)}=0.011010011001{\dots }_{(2)}$

    $⟹1678.4125_{(10)}=11010001110.0110100110011001{\dots }_{(2)}$

（2）将$11010001110.0110100110011001\dots$改写成2进制科学记数法：
$$11010001110.0110100110011001\dots =1.10100011100110100110011001\dots \times 2^{10}$$

（3）将 $1.10100011100110100110011001\dots \times 2^{10}$ 存入64位二进制码中。

  说明：

  为达到最大限度利用64位二进制码，我们只存储最有效的部分：
  ①：符号值 sign
  ②：指数值 Exponent（显然不用存储底数2，因为二进制必定是2）
  ③：小数值 Mantissa（显然不用存储整数部分1，因为科学技术法必定有整数部分1）

现在来规划32位码，将其分割成3个部分以存储上述3个重要内容，如下：

  ① 符号位：
  1 bit：0表示正数，1表示负数

  ② 指数位：
  8 bit：取值 [ 0，255 ] 。但是为了表示负指数，所以以 127 作为计算分界点，将取值范围变为 [ -127，127 ] 。也就是说：在二进制转浮点数时，需要将二进制码的指数值减 127 才能得到真实的指数值。如：$1.2\times 2^4$，在二进制指数位实际存储的是127+4的二进制码，也就是10000011。

  ③ 小数位:
  23 bit： 直接放入二进制小数位即可

  实际存入：
$$+1.10100011100110100110011001\dots \times 2^{10}$$

注意：

  这里就产生了截断误差：

  原本小数位是：$10100011100110100110011001\dots$

  现在小数位是：$10100011100110100110011$

  如果转换回来：

  二进制码：     $01000100110100011100110100110011$

  二进制科学计数法： $1.10100011100110100110011\times 2^{10}$

  二进制：      $11010001110.0110100110011$

  十进制：      $1678.4124755859375$

  差值:        $|1678.4125-1678.4124755859375|=0.0000244140625=2.44\times 10^{-5}<10^{-4}$

所以在浮点数判断时，通常写作:

//仅针对 32 位架构，当代 64 位架构缩小到 e-6
if( abs(a,b) < e-4 )
{
	printf("float a is equal to float b\n");
}

---

（二） IEEE 754 浮点数表示法（二进制码转数值）

在第一节中，我们了解了IEEE 754二进制码由
$$SignExponentMantissa$$

组成，那么计算机如何将这3个二进制数 $S$、$E$、$M$ 还原成数值 $V$？

1. 以M和E代表数值

根据 IEEE 754 原理，我们逆向可得公式：

$$V=(1+M\times 2^{23})\times 2^{E-127}$$

即：

$$V=(1+\frac{M}{2^{23}})\times 2^{E-127}$$

我们来验证一下，还是以$1678.4125_{(10)}$ 举例：

10 进制数值：$1678.4125$

2 进制数值：$11010001110.0110100110011001\dots$

2 进制码：$01000100110100011100110100110011$

则，我们可以得到3个二进制数：

① $S=0$

② $E=10001001$

③ $M=10100011100110100110011$

$\because M=10100011100110100110011$

$⟹\frac{M}{2^{23}}=0.10100011100110100110011$

$⟹1+\frac{M}{2^{23}}=1.10100011100110100110011$

$\because E=10001001_{(2)}$

$⟹E-127_{(10)}=137-127=10$

$⟹2^{E-127_{(10)}}=2^{10}$

$\therefore V=(1+\frac{M}{2^{23}})\times 2^{E-127}=1.10100011100110100110011\times 2^{10}$

$\therefore V=11010001110.0110100110011_{(2)}$

$\therefore V=1678.4124755859375_{(10)}$

（在第一节说明了这里产生了截断误差）

2. 以M和E代表二进制码

$$V=M+2^{23}E$$

还是以$1678.4125_{(10)}$ 举例：

10 进制数值：$1678.4125$

2 进制数值：$11010001110.0110100110011001\dots$

2 进制码：$01000100110100011100110100110011$

则，我们可以得到3个二进制数：

① $S=0$

② $E=10001001$

③ $M=10100011100110100110011$

那么：

$\begin{array}{rlr}{2}^{23}E& =10001001& 00000000000000000000000\\ M& =& 10100011100110100110011\\ M+2^{23}E& =10001001& 10100011100110100110011\end{array}$

$$⟹V=M+2^{23}E$$

3. 重要推论：

$$V=M+2^{23}E=(1+\frac{M}{2^{23}})\times 2^{E-127}$$

---

（三） IEEE 754 浮点数表示法（神奇的性质）

既然:

$$V=(1+\frac{M}{2^{23}})\times 2^{E-127}$$

那么对等式左右两边取对数可得：

$$\begin{array}{rl}{\log }_{2}V& ={\log }_2[(1+\frac{M}{2^{23}})\times 2^{E-127}]\\ & ={\log }_2[(1+\frac{M}{2^{23}})]+{\log }_2{2}^{E-127}]\\ & ={\log }_2[(1+\frac{M}{2^{23}})]+E-127\end{array}$$

$\because M\in [0,1.1111111111111111111111\times 2^{22}]$

$⟹\frac{M}{2^{23}}\in [0,1)$

$\because {\log }_2(1+x)\approx x+\delta$

$\begin{array}{rl}\therefore {\log }_2[(1+\frac{M}{2^{23}})]+E-127& \approx \frac{M}{2^{23}}+\delta +E-127\\ & \approx \frac{M}{2^{23}}+E-127+\delta \\ & \approx \frac{1}{2^{23}}(M+2^{23}E)-127+\delta \end{array}$

所以得到一个重要推论：

$$lo{g}_{2}V=\frac{1}{2^{23}}V-127+\delta$$

---

（四）求解

要求 $y=\frac{1}{\sqrt{x}}$

即求 $y=x^{-\frac{1}{2}}$