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

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

简介

介绍使用魔数0x2a517d47快速求立方根$\sqrt[3]{x}$的C语言实现和公式的推导。

代码

float MagicCubeRoot(float x)
{
    float xthird = 0.333f * x;
    int i = *(int*)&x;
    i = (0x2a517d47 + (0.333f * i));
    x = *(float*)&i;
    x = 0.667f * x + xthird / (x * x);
    
    return x;
}

代码用于快速计算立方根$\sqrt[3]{x}$。

代码中核心部分是

i = (0x2a517d47 + (0.333f * i));

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

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

x = 0.667f * x + xthird / (x * x);

所以整个计算过程就是

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

将输入的数转换成整数

2.i = (0x2a517d47 + (0.333f * i));

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

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

转换回浮点数。

4.x = 0.667f * x + xthird / (x * x);

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

完成快速计算立方根$\sqrt[3]{x}$。

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

float MagicCubeRoot(float x)
{
    float xthird = 0.333f * x;
    int i = *(int*)&x;
    i = (0x2a517d47 + (0.333f * i));
    x = *(float*)&i;
    x = 0.667f * x + xthird / (x * x);
    x = 0.667f * x + xthird / (x * x);
    
    return x;
}

注意，上面的函数只支持正数，传进来的值x需要大于或等于0。

如果要支持负数，可以增加如下判断处理。

float MagicCubeRoot(float x)
{
    if (x < 0)
    {
        x = -x;
        
        float xthird = 0.333f * x;
        int i = *(int*)&x;
        i = (0x2a517d47 + (0.333f * i));//i = (int) (0x2a517d3c + (0.333f * i));
        x = *(float*)&i;
        x = 0.667f * x + xthird / (x * x);
        //x = 0.667f * x + xthird / (x * x);
        
        x = -x;
    }
    else
    {
        float xthird = 0.333f * x;
        int i = *(int*)&x;
        i = (0x2a517d47 + (0.333f * i));//i = (int) (0x2a517d3c + (0.333f * i));
        x = *(float*)&i;
        x = 0.667f * x + xthird / (x * x);
        //x = 0.667f * x + xthird / (x * x);
    }
    
    return x;
}

牛顿迭代法计算方式

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

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

我们计算立方根的公式：

$y=\sqrt[3]{x}=x^{\frac{1}{3}}$

所以

$y^3=x$

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

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

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

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

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

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

所以

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

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

$y_2=(2\times 2+2÷2^2)÷3=1.5$

$y_3=(2\times 1.5+2÷1.5^2)÷3=1.296296$

$y_4=(2\times 1.296296+2÷1.296296^2)÷3=1.260932$

$y_5=(2\times 1.260932+2÷1.260932^2)÷3=1.259922$

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

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

算法讲解

接下来解释该算法的核心i = (0x2a517d47 + (0.333f * i));。

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

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[3]{x}=x^{\frac{1}{3}}$

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

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

将x和y用浮点数替换：

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

${\log }_2(1+m_y)+e_y=\frac{1}{3}({\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}{3}({\log }_2(1+m_x)+e_x)$

进行化简

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

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

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

进行化简

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

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

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

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

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

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

用代码表示就是：

i = K + i / 3;

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

魔数的选定

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

关于$\sigma$值的选定的详细内容，可以参考以下论文：

Lomont, Chris. Fast Inverse Square Root. February 2003.

McEniry, Charles. The Mathematics Behind the Fast Inverse Square Root Function Code . August 2007 .

 

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

带入计算出K值：

$K=\frac{2}{3}L(B-\sigma )=\frac{2}{3}{2}^{23}(127-0.0450465)\approx 709983559=0x2a517d47$

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

---

参考资料

平方根倒数速算法
0x5f3759df

---

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