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

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

简介

介绍使用魔数0x1fbd1df5快速求平方根$\sqrt{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;
}

代码用于快速计算平方根$\sqrt{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;

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

完成快速计算平方根$\sqrt{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;
    }
}

牛顿迭代法计算方式

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

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

我们计算平方根的公式：

$y=\sqrt{x}=x^{\frac{1}{2}}$

所以

$y^2=x$

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

$f(y)=y^2-x=0$

$f^{\prime }(y)=2y=0$

将$f(y)$和$f^{\prime }(y)$带入

$y-\frac{f(y)}{f^{\prime }(y)}=y-\frac{y^2-x}{2y}$

​ $=\frac{1}{2}(y+\frac{x}{y})$

所以

$y_{n+1}=\frac{1}{2}(y_n+\frac{x}{y_n})$ = $0.5y_n+0.5x÷y_n$

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

$y_2=0.5\times 2+0.5\times 2÷2=1.5$

$y_3=0.5\times 1.5+0.5\times 2÷1.5=1.416667$

$y_4=0.5\times 1.416667+0.5\times 2÷1.416667=1.414216$

$y_5=0.5\times 414216+0.5\times 2÷414216=1.414214$

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

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

算法讲解

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

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

32位浮点数基本概念

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

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

用公式表示就是

$(-1{)}^s(1+m)2^e$

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

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

$(1+m)2^e$

指数部分$e$的数值范围，$-127\le e\le 128$。

尾数部分$m$的数值范围，$0\le m<1$。

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

$M+LE$

这里$E$表示指数(exponent)，$M$表示尾数(mantissa)，$L$是$2^{23}$。

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

所以$E$和$e$的转换关系是$e=E-B$，$B$是127。

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

所以$M$和$m$的转换关系是$m=\frac{M}{L}$，$L$是$2^{23}$。

推导过程

接下来进入推导过程。

给定一个$x$，计算平方根$y$：

$y=\sqrt{x}=x^{\frac{1}{2}}$

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

${\log }_{2}y=\frac{1}{2}{\log }_{2}x$

将x和y用浮点数替换：

${\log }_2((1+m_y)2^{e_y})=\frac{1}{2}({\log }_2((1+m_x)2^{e_x}))$

${\log }_2(1+m_y)+e_y=\frac{1}{2}({\log }_2(1+m_x)+e_x)$

算式两边都有这样的项

${\log }_2(1+v)$

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

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

方程式：

${\log }_2(1+v)\approx v+\sigma$

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

上面的方程

${\log }_2(1+m_y)+e_y=\frac{1}{2}({\log }_2(1+m_x)+e_x)$

进行化简

$m_y+\sigma +e_y\approx \frac{1}{2}(m_x+\sigma +e_x)$

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

$\frac{M_y}{L}+\sigma +E_y-B\approx \frac{1}{2}(\frac{M_x}{L}+\sigma +E_x-B)$

进行化简

$\frac{M_y}{L}+E_y\approx \frac{1}{2}(\frac{M_x}{L}+\sigma +E_x-B)-\sigma +B$

$\frac{M_y}{L}+E_y\approx \frac{1}{2}(\frac{M_x}{L}+E_x)-\frac{1}{2}(\sigma -B)$

$M_y+L{E}_y\approx \frac{1}{2}L(B-\sigma )+\frac{1}{2}(M_x+L{E}_x)$

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

${\mathbf{I}}_{\mathbf{y}}\approx \frac{1}{2}L(B-\sigma )+\frac{1}{2}{\mathbf{I}}_{\mathbf{x}}$

其中$\frac{1}{2}L(B-\sigma )$是一个常数。

用代码表示就是：

i = K + (i >> 1);

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

魔数的选定

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

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

带入计算出K值：

$K=\frac{1}{2}L(B-\sigma )=\frac{1}{2}{2}^{23}(127-0.0450465)\approx 532487669=0x1fbd1df5$

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

参考资料

平方根倒数速算法
0x5f3759df

---

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