对于快速计算出平方根倒数 $\frac{1}{\sqrt{x}}$，有一个上世纪90年代的经典算法：

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;
}

扫盲：浮点数在计算机中的储存

float32和double64都由 IEE754 标准定义，在这里只简单了解float32。

在这32位中分为三个部分:

• Sign 符号：0代表整数，1代表负数
• Exponent 指数位

二进制的8bit可以表示256种状态，而IEE754规定，这八位用于表示[-127, 127]范围内的指数。 为了方便地表示负数，其规定在float32中指数的偏移量为127，因此指数的实际值为e = exponent - 127。 这将确保以二进制储存的biased exponent总是非负整数。

• Fraction 尾数位

尾数位在格式中有23位，因此可以表示的精度为2^-23，约等于1.19 * 10^-7。

如有一个十进制数13.62，将整数部分除2取余逆序排列为1101，将小数部分乘2取整顺序排列为101。这个十进制小数在二进制的表示为1101.101，即1.101101 * 2^3。

将指数部分加上偏移量127，得到130，则exponent部分的值为10000010。

因为规定了小数的最高位总是非零数，所以在二进制中它总是为1，因此在尾数部分可以省略最高位的储存，只考虑小数点后面的数字。如凑不够23位则低位补零。

$$ {\text{value}}=(-1)^{\text{sign}}\times 2^{(E-127)}\times \left(1+\frac{M}{2^{23}}\right) $$

其中，$\frac{M}{2^{23}}$ 将尾数 $M$ 归一化到范围 $[0, 1)$。

牛顿迭代法求根

对于 $ f(x)=0 $ 有一个近似解 $ x_n $，通过给出的根得到更近似的根。

$$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} $$

已经证明牛顿迭代法的二次收敛条件： ${\displaystyle f'(x)\neq 0}$ ; ${\displaystyle x\in I}$，其中 ${\displaystyle I}$ 为区间$[α − r, α + r]$ ; ${\displaystyle f(x)}$ 在 ${\displaystyle I}$ 上连续可微; 初始根 $x_0$ 足够接近根 $x_*$

最开始的浮点数公式： $x=\left(1+\frac{M}{2^{23}}\right) 2^{(E-127)}$

两边同取对数，

$$ \begin{align} \log_{2}(x) &= \log_{2}(\left(1+\frac{M}{2^{23}}\right) 2^{(E-127)}) \\ &= \left(E-127\right)+\log_{2}(1+\frac{M}{2^{23}}) \end{align} $$

在$ x \in [0, 1]$，$\log_{2}(1+x) \approx x$，$ y=\log_{2}(1+x) $ 与 $y=x$ 在图像上十分接近。

$$ \begin{align} \log_{2}(x) &= \left(E-127\right)+\log_{2}(1+\frac{M}{2^{23}})\\ &= \frac{M}{2^{23}}+E-127 \\ &= \frac{M+2^{23} \times E}{2^{23}}-127 \end{align} $$

此时，$M+2^{23}E$ 就是浮点数在二进制中的表达，其中 $2^{23}$ 使 $E$ 在二进制中向左移23位。

$log_{2}(x) = \frac{M+2^{23} \times E}{2^{23}}-127 $ 其实就是浮点数 $x$ 与其在二进制中的关系

设有解 $a=\frac{1}{\sqrt{y}}$

$\log_{2}(a) = \log_{2}(\frac{1}{\sqrt{y}}) = -\frac{1}{2}\log_{2}(y)$

将浮点数 $a$ 和 $y$ 转化为二进制形式 $A, Y$，代入上式，得到

$\frac{A}{2^{23}}-127 = -\frac{1}{2}\left( \frac{Y}{2^{23}}-127 \right) $

$A = 381 \times 2^{22} + α - \frac{1}{2}Y$

这个式子，其实就是这行代码的含义：

i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 

但 0x5f3759df 这个值是怎么来的呢？或句话说，如何计算出最佳的修正因子$α$？

引入一个修正因子$α$，使得直线 $y=x$ 上移与 $y=\log_{2}(1+x)$ 在图像上更加接近。

$$\log_{2}(1+x) \approx x+α$$

切比雪夫最佳逼近

$$ E(\alpha) = \max_{x \in [0,b]} \left| (x + \alpha) - \log_2(1 + x) \right| $$

最佳逼近直线为

$$y=x+0.0431$$

Refs

float32 Picture (By Fresheneesz at the English Wikipedia project, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=3357169)

NewtonMethod Picture 作者 Ralf Pfeifer - de:Image:NewtonIteration Ani.gif，CC BY-SA 3.0，https://commons.wikimedia.org/w/index.php?curid=2268473