当前位置 博文首页 > dadalaohua的博客:【学习笔记】使用魔数快速求平方根

    dadalaohua的博客:【学习笔记】使用魔数快速求平方根

    作者:[db:作者] 时间:2021-07-27 14:44

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

    简介

    介绍使用魔数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;
        }
    }
    

    牛顿迭代法计算方式

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

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

    我们计算平方根的公式:

    y = x = x 1 2 y = {\sqrt{x}} = x^{\frac 12} y=x ?=x21?

    所以

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

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

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

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

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

    y ? f ( y ) f ′ ( y ) = y ? y 2 ? x 2 y y - \frac{f(y)}{f'(y)} = y - \frac{{{y}^{2}} - x}{{2}{y}}