← 返回全部文章

快速求解平方根倒数算法

[ Code & Dev ]
[*] 由AI生成的摘要 GPT

本文介绍了一种快速计算平方根倒数的算法，该算法源于上世纪90年代。文章首先解释了浮点数在计算机中的存储方式，特别是float32格式的结构，包括符号位、指数位和尾数位。接着，介绍了牛顿迭代法的基本原理及其在求解平方根倒数中的应用。通过对浮点数的对数变换，推导出与平方根倒数相关的公式，并解释了代码中使用的神秘常数 0x5f3759df 的来源，最后提到切比雪夫最佳逼近的概念，以优化计算结果。

对于快速计算出平方根倒数 1x\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。

float

在这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位则低位补零。

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

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

牛顿迭代法求根

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

xn+1=xn−f(xn)f′(xn)x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

NewtonMethod

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

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

两边同取对数，

log⁡2(x)=log⁡2((1+M223)2(E−127))=(E−127)+log⁡2(1+M223)\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∈[0,1] x \in [0, 1]，log⁡2(1+x)≈x\log_{2}(1+x) \approx x，y=log⁡2(1+x)y=\log_{2}(1+x) 与 y=xy=x 在图像上十分接近。

log⁡2(x)=(E−127)+log⁡2(1+M223)=M223+E−127=M+223×E223−127\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+223EM+2^{23}E 就是浮点数在二进制中的表达，其中 2232^{23} 使 EE 在二进制中向左移23位。

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

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

log⁡2(a)=log⁡2(1y)=−12log⁡2(y)\log_{2}(a) = \log_{2}(\frac{1}{\sqrt{y}}) = -\frac{1}{2}\log_{2}(y)

将浮点数 aa 和 yy 转化为二进制形式 A,YA, Y，代入上式，得到

A223−127=−12(Y223−127)\frac{A}{2^{23}}-127 = -\frac{1}{2}\left( \frac{Y}{2^{23}}-127 \right)

A=381×222+α−12YA = 381 \times 2^{22} + α - \frac{1}{2}Y

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

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

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

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

log⁡2(1+x)≈x+α\log_{2}(1+x) \approx x+α

切比雪夫最佳逼近

E(α)=max⁡x∈[0,b]∣(x+α)−log⁡2(1+x)∣E(\alpha) = \max_{x \in [0,b]} \left| (x + \alpha) - \log_2(1 + x) \right|

最佳逼近直线为

y=x+0.0431y=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