快速计算平方根倒数

by · 2019年11月27日 · 2382 Words · ~5min reading time | Improve on

【学习笔记】快速平方根倒数算法 - w450468524的专栏 - CSDN博客

很以前久，看《DOOM 启示录》的时候就看到，当年卡马克大神在《雷神之锤》中使用了一个神奇的数字，能够通过位操作快速计算平方根倒数 y=1x√y=1xy=\frac{1}{\sqrt{x}}。但是当时并没有深究，这两天偶然看到了这篇文章，终于将我这个多年的 “未解之谜” 解开了，特此做下笔记，图片取自原文。

代码实现

首先是这个神奇算法的代码：

float fast_inverse_sqrt(float x)
{
    float half_x = 0.5 * x;
    int i = *((int *)&x); // 以整数方式读取X
    i = 0x5f3759df - (i>>1); // 神奇的步骤
    x = *((float *)&i); // 再以浮点方式读取i
    x = x*(1.5 - (half_x * x * x)); // 牛顿迭代一遍提高精度
    return x;
}

其中最最重要的就是第二步，通过这个神奇的步骤，能够得到 x 平方根倒数的近似值。前后穿插着用整型读取浮点数等等，实在是不知所云。但实际上，其中包含着很有意思的数学背景（每次一提到数学，我就觉得好方）。

背后的数学

首先来重新温习一下机组里面的浮点数基础知识（机组 60 分默默飘过）。

单精度浮点数

单精度浮点数有 32bit，包括 1 位符号位 s，8 位指数位 e，23 位尾数位 m（赞一下 CSDN 的 markdown 编辑器的 latex 数学公式的功能~)：

s e⋯e⏟8 m⋯m⏟23se⋯e⏟8m⋯m⏟23\qquad s \space \underbrace{e\cdots e}8\space \underbrace{m\cdots m}{23}

其表示的数值为：(1+m′)2e′(1+m′)2e′(1+m')2^{e'}，这里的 m’代表尾数部分的数值，m’只包含小数位，即 m′=0.m⋯m⏟23m′=0.m⋯m⏟23m'=0.\underbrace{m\cdots m}_{23 }，e’则是指数部分的数值，只有整数位，e′=e⋯e⏟8−Be′=e⋯e⏟8−Be'=\underbrace{e\cdots e}_8-B。为了表示负的指数，指数部分要减去一个偏移量 B。
再规定 M 表示以整数方式解释尾数二进制位的值，E 表示相同的方式解释指数部分的值。可以得到以下的转换关系：

m′=MLL=223m′=MLL=223\qquad m'=\frac{M}{L}\qquad L=2^{23}
e′=E−BB=127e′=E−BB=127\qquad e'=E-B\qquad B=127

那么对于一个浮点数的二进制位 s e⋯e⏟8 m⋯m⏟23se⋯e⏟8m⋯m⏟23s \space \underbrace{e\cdots e}8\space \underbrace{m\cdots m}{23}以整数方式解读的话，其值就为：I=EL+MI=EL+MI=EL+M（因为平方根输入只能为正，所以默认 s 为 0）。

平方根倒数近似计算

下面开始进入正题：
对于输入 x，我们要计算的是 y=1x√y=1xy=\frac{1}{\sqrt{x}}，将该式子两边取对数：log2(y)=−12log2(x)log2⁡(y)=−12log2⁡(x)\log_2(y)=-\frac{1}{2}\log_2(x)。然后我们将 x，y 都替换为浮点数表示形式：

log2((1+m′y)2e′y)=−12log2((1+m′x)2e′x)log2⁡((1+my′)2ey′)=−12log2⁡((1+mx′)2ex′)\qquad \log_2((1+m'_y)2^{e'_y})=-\frac{1}{2}\log_2((1+m'_x)2^{e'_x})
log2(1+m′y)+e′y=−12(log2(1+m′x)+e′x)log2⁡(1+my′)+ey′=−12(log2⁡(1+mx′)+ex′)\qquad \log_2(1+m'_y)+e'_y=-\frac{1}{2}(\log_2(1+m'_x)+{e'_x})

观察发现，上式等号两边均有 log2(1+v)log2⁡(1+v)\log_2(1+v)形式的式子，v 的范围是 0 到 1，而该式在这个定义域内的函数图像非常接近一条直线：

因此我们近似的用 x+σx+σx+\sigma替换得到：

m′y+σ+e′y=−12(m′x+σ+e′x)my′+σ+ey′=−12(mx′+σ+ex′)\qquad m'_y+\sigma+e'_y=-\frac{1}{2}(m'_x+\sigma+{e'_x})

σσ\sigma为预先计算的一个常数。接下来将数值转为整形方式读取二进制位的形式：

MyL+σ+Ey−B=−12(MxL+σ+Ex−B)MyL+σ+Ey−B=−12(MxL+σ+Ex−B)\qquad \frac{M_y}{L}+\sigma+E_y-B=-\frac{1}{2}(\frac{M_x}{L}+\sigma+E_x-B)

移项合并后得到：

32L(σ−B)+My+LEy=−12(Mx+LEx)32L(σ−B)+My+LEy=−12(Mx+LEx)\qquad \frac{3}{2}L(\sigma-B)+M_y+LE_y=-\frac{1}{2}(M_x+LE_x)

恰好我们又发现，在等式的两边都含有以整数方式解释浮点数的表示 I=M+LEI=M+LEI=M+LE，进一步化简有：

Iy=−12Ix+κκ=32L(B−σ)Iy=−12Ix+κκ=32L(B−σ)\qquad I_y=-\frac{1}{2}I_x + \kappa\qquad \kappa=\frac{3}{2}L(B-\sigma)

由此，我们得到了一个比较有用的近似公式，κκ\kappa便是神奇的数字 0x5f3759df 。该公式表明计算 x 的平方根倒数时，我们只需要用整型方式读取 x，再代入上式计算得到结果的整形解释值，然后再用浮点数方式读取即可！我当时就震惊了。

牛顿迭代提高精度

经过上一步，我们有了初步的近似结果，为了进一步提高精度，我们再用牛顿迭代法计算一次。牛顿迭代法描述的是，给定连续函数 f(x)，对于方程 f(x)=0，确定任意初始的 x0x0x_0，我们可以通过下面这个迭代式计算得到其解：

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

我们要计算的式子为 y=1x0√，x0y=1x0，x0y=\frac{1}{\sqrt{x_0}}，x_0为已知（这个 x0x0x_0是指输入，不是牛顿迭代的初始值），将其转换为关于 y 的方程就是 f(y)=1y2−x0=0f(y)=1y2−x0=0f(y)=\frac{1}{y^2}-x_0=0，根据牛顿迭代可得：

yn+1=yn−1y2n−x0−2y3nyn+1=yn−1yn2−x0−2yn3\qquad y_{n+1}=y_n-\frac{\frac{1}{y_n^2}-x_0}{\frac{-2}{yn^3}}
yn+1=32yn−12x0y3nyn+1=32yn−12x0yn3\qquad y{n+1}=\frac{3}{2}y_n-\frac{1}{2}x_0y_n^3

也就是最后一步的代码了。

快速计算平方根倒数

代码实现

背后的数学

单精度浮点数

平方根倒数近似计算

牛顿迭代提高精度

更多讨论

分类

快速计算平方根倒数

代码实现

背后的数学

单精度浮点数

平方根倒数近似计算

牛顿迭代提高精度

更多讨论

分类

标签云