【学习笔记】使用魔数快速求平方根

简介

介绍使用魔数0x1fbd1df5快速求平方根x{\sqrt{x}}x的C语言实现和公式的推导。

代码

float MagicSqrt(float x)
{float xhalf = 0.5f * x;int i = *(int*)&x;i = 0x1fbd1df5 + (i >> 1);x = *(float*)&i;x = 0.5f * x + xhalf / x;return x;
}

代码用于快速计算平方根x{\sqrt{x}}x。

代码中核心部分是

i = 0x1fbd1df5 + (i >> 1);

该行代码就完成了计算平方根。

另外再使用一次牛顿迭代法提高下精度

x = 0.5f * x + xhalf / x;

所以整个计算过程就是

1.i = *(int*)&x;

将输入的数转换成整数

2.i = 0x1fbd1df5 + (i >> 1);

通过魔数完成平方根的计算。

3.x = *(float*)&i;

转换回浮点数。

4.x = 0.5f * x + xhalf / x;

使用一次牛顿迭代法提高下精度。

完成快速计算平方根x{\sqrt{x}}x。

如果需要提高精度，可以多进行一次牛顿迭代。

float MagicSqrt(float x)
{float xhalf = 0.5f * x;int i = *(int*)&x;i = 0x1fbd1df5 + (i >> 1);x = *(float*)&i;x = 0.5f * x + xhalf / x;x = 0.5f * x + xhalf / x;return x;
}

传进来的值x必须大于或等于0。

可以加入一个判断防止传入的值小于0.

float MagicSqrt(float x)
{if (x < 0){return -1;}else{float xhalf = 0.5f * x;int i = *(int*)&x;i = 0x1fbd1df5 + (i >> 1);x = *(float*)&i;x = 0.5f * x + xhalf / x;//x = 0.5f * x + xhalf / x;return x;}
}

牛顿迭代法计算方式

一般计算平方根可以使用牛顿迭代法。

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

我们计算平方根的公式：

y=x=x12y = {\sqrt{x}} = x^{\frac 12}y=x=x21

所以

y2=x{y}^2 = xy2=x

构建以y为自变量的函数方程为

f(y)=y2−x=0f(y) = {y}^2 - x = 0f(y)=y2−x=0

f′(y)=2y=0f'(y) = {2}{y} = 0f′(y)=2y=0

将f(y)f(y)f(y)和f′(y)f'(y)f′(y)带入

y−f(y)f′(y)=y−y2−x2yy - \frac{f(y)}{f'(y)} = y - \frac{{{y}^{2}} - x}{{2}{y}}y−f′(y)f(y)=y−2yy2−x

=12(y+xy)= \frac{1}{2}{(y + \frac{x}{y})}=21(y+yx)

所以

yn+1=12(yn+xyn)y_{n+1} = \frac{1}{2}{(y_{n} + \frac{x}{y_{n}})}yn+1=21(yn+ynx) = 0.5yn+0.5x÷yn0.5{y_n} + 0.5{x}÷{y_n}0.5yn+0.5x÷yn

实际计算，设x=2x = 2x=2，y1=2y_1 = 2y1=2

y2=0.5×2+0.5×2÷2=1.5y_2 = 0.5 × 2 + 0.5 × 2 ÷ 2 = 1.5y2=0.5×2+0.5×2÷2=1.5

y3=0.5×1.5+0.5×2÷1.5=1.416667y_3 = 0.5 × 1.5 + 0.5 × 2 ÷ 1.5 = 1.416667y3=0.5×1.5+0.5×2÷1.5=1.416667

y4=0.5×1.416667+0.5×2÷1.416667=1.414216y_4 = 0.5 × 1.416667 + 0.5 × 2 ÷ 1.416667 = 1.414216y4=0.5×1.416667+0.5×2÷1.416667=1.414216

y5=0.5×414216+0.5×2÷414216=1.414214y_5 = 0.5 × 414216 + 0.5 × 2 ÷ 414216 = 1.414214y5=0.5×414216+0.5×2÷414216=1.414214

经过多次迭代，计算出来的结果2≈1.414214{\sqrt{2}} \approx 1.4142142≈1.414214。

使用牛顿迭代法计算平方根需要大量的浮点数运算，运算速度会远远弱于整数运算。

算法讲解

接下来解释该算法的核心i = 0x1fbd1df5 + (i >> 1);。

首先需要知道浮点数的基本概念。

32位浮点数基本概念

浮点数由三部分组成：符号，指数和尾数。

32位浮点数，用二进制表示

用公式表示就是

(−1)s(1+m)2e(-1)^s(1+m)2^e(−1)s(1+m)2e

这里s是符号(sign)，e是指数(exponent)，m是尾数(mantissa)。

因为这里计算平方根只对正数有意义，所以从假设符号位是 0 ，公式简化为

(1+m)2e(1+m)2^e(1+m)2e

指数部分eee的数值范围，−127≤e≤128-127 \leq e \leq 128−127≤e≤128。

尾数部分mmm的数值范围，0≤m<10 \leq m < 10≤m<1。

如果将浮点数转换成整数时，整数的数值就是

M+LEM + LEM+LE

这里EEE表示指数(exponent)，MMM表示尾数(mantissa)，LLL是2232^{23}223。

将指数部分看做是整数，用EEE来表示，那么范围是0≤E≤2550 \leq E \leq 2550≤E≤255。EEE如果减去127，范围变成-127 - 128，变成一个有符号数。

所以EEE和eee的转换关系是e=E−Be = E - Be=E−B，BBB是127。

如果将尾数部分看做是整数，用MMM来表示。数值范围是0≤M<2230 \leq M < 2^{23}0≤M<223。

所以MMM和mmm的转换关系是m=MLm = \frac{M}{L}m=LM，LLL是2232^{23}223。

推导过程

接下来进入推导过程。

给定一个xxx，计算平方根yyy：

y=x=x12y = {\sqrt{x}} = x^{\frac 12}y=x=x21

首先对等式两边取以 2 为底的对数

log⁡2y=12log⁡2x\log_2 y = {\frac 12}\log_2 xlog2y=21log2x

将x和y用浮点数替换：

log⁡2((1+my)2ey)=12(log⁡2((1+mx)2ex))\log_2 ((1+m_y)2^{e_y}) = {\frac 12}(\log_2 ((1+m_x)2^{e_x}))log2((1+my)2ey)=21(log2((1+mx)2ex))

log⁡2(1+my)+ey=12(log⁡2(1+mx)+ex)\log_2 (1+m_y) + e_y = {\frac 12}(\log_2 (1+m_x) + e_x)log2(1+my)+ey=21(log2(1+mx)+ex)

算式两边都有这样的项

log⁡2(1+v)\log_2(1 + v)log2(1+v)

其中v的值范围0 < vvv < 1。

当v的取值在0到1之间时，这个函数和一条直线很接近：

方程式：

log⁡2(1+v)≈v+σ\log_2(1 + v) \approx v + \sigmalog2(1+v)≈v+σ

其中σ\sigmaσ是一个常数，可以通过调整这个常数让两个曲线更加近似。

上面的方程

log⁡2(1+my)+ey=12(log⁡2(1+mx)+ex)\log_2 (1+m_y) + e_y = {\frac 12}(\log_2 (1+m_x) + e_x)log2(1+my)+ey=21(log2(1+mx)+ex)

进行化简

my+σ+ey≈12(mx+σ+ex)m_y + \sigma + e_y \approx {\frac 12}(m_x + \sigma + e_x)my+σ+ey≈21(mx+σ+ex)

接下来用整形解释下的指数和尾数来替代浮点数解释：

MyL+σ+Ey−B≈12(MxL+σ+Ex−B)\frac{M_y}{L} + \sigma + E_y - B \approx {\frac 12}(\frac{M_x}{L} + \sigma + E_x - B)LMy+σ+Ey−B≈21(LMx+σ+Ex−B)

进行化简

MyL+Ey≈12(MxL+σ+Ex−B)−σ+B\frac{M_y}{L} + E_y \approx {\frac 12}(\frac{M_x}{L} + \sigma + E_x - B) - \sigma + BLMy+Ey≈21(LMx+σ+Ex−B)−σ+B

MyL+Ey≈12(MxL+Ex)−12(σ−B)\frac{M_y}{L} + E_y \approx {\frac 12}(\frac{M_x}{L} + E_x) - \frac{1}{2}(\sigma - B)LMy+Ey≈21(LMx+Ex)−21(σ−B)

My+LEy≈12L(B−σ)+12(Mx+LEx)M_y + LE_y \approx {\frac 12}L(B - \sigma) + {\frac 12}(M_x + LE_x)My+LEy≈21L(B−σ)+21(Mx+LEx)

两边都得到了整形解释的值：

Iy≈12L(B−σ)+12Ix\mathbf{I_y} \approx {\frac 12}L(B - \sigma) + {\frac 12}\mathbf{I_x}Iy≈21L(B−σ)+21Ix

其中12L(B−σ){\frac 12}L(B - \sigma)21L(B−σ)是一个常数。

用代码表示就是：

i = K + (i >> 1);

K就是算法中的常数。通过选取合适的σ\sigmaσ的值，就能得到K值。

魔数的选定

计算魔数所用的可以采用穷举搜索或者二分法寻找最优值。

这里选定的值为σ\sigmaσ = 0.0450465。

带入计算出K值：

K=12L(B−σ)=12223(127−0.0450465)≈532487669=0x1fbd1df5K = {\frac 12}L(B - \sigma) = {\frac 12}{2^{23}}(127 - 0.0450465) \approx 532487669 = 0x1fbd1df5K=21L(B−σ)=21223(127−0.0450465)≈532487669=0x1fbd1df5

这里的K值就是公式中的0x1fbd1df5。

参考资料

平方根倒数速算法
0x5f3759df

本文链接：https://blog.csdn.net/u012028275/article/details/113793827