不用sqrt()函数，求平方根的三种方法

最新推荐文章于 2023-03-31 00:26:37 发布

Donald Su

最新推荐文章于 2023-03-31 00:26:37 发布

阅读量4.7k

点赞数 3

分类专栏： C/C++ LeetCode 文章标签：面试

本文链接：https://blog.csdn.net/donaldsy/article/details/106339865

版权

LeetCode 同时被 2 个专栏收录

30 篇文章 0 订阅

订阅专栏

C/C++

11 篇文章 0 订阅

订阅专栏

最近看到了这个比较有意思的题目，探究了一下。

文章目录

当然有最暴力的方法，直接遍历，(0.0, x)区间内所有的数据(也可以是 x/2)，看值是否相等，但该方法太过暴力，在此不讨论。

可以采用下面的 二分法、牛顿法，Quake III源码中的方法。

浮点数的比较方法比较特别，不能使用平常的 x == 0的方法进行比较，下面都用到了equal比较方法，equal的代码：

const double Min_value = 1e-7; //定义最小的浮点数数误差，这里的用法比 宏更好
bool equal(double num1, double num2)
{
    if ((num1 - num2 > -Min_value) && (num1 - num2 < Min_value))
        return true;
    else
        return false;
}

1、二分法

类似于二分查找的思路，时间复杂度 $O (l o g n)$ 。

double sqrt_binary_search(double square)
{
    if(square < Min_value)
        return -1;

    double low = 0.0;
    double high = square;

    // while (high - low > Min_value)
    while (!equal(high, low))
    {
        /* code */
        double mid = (high + low) / 2;
        if (mid * mid > square)
        {
            high = mid;
        }
        else
        {
            low = mid;
        }
    }

    return low;
}

2、牛顿法

先分析下牛顿法的原理，用的是牛顿迭代方法。

牛顿迭代法又称为牛顿-拉弗森方法，实际上是由牛顿、拉弗森各自独立提出来的。牛顿-拉弗森方法提出来的思路就是利用切线是曲线的线性逼近这个思想。

随便找一个曲线上的A点（为什么随便找，根据切线是切点附近的曲线的近似，应该在根点附近找，但是很显然我们现在还不知道根点在哪里），做一个切线，切线的根（就是和x轴的交点）与曲线的根，还有一定的距离。牛顿、拉弗森们想，没关系，我们从这个切线的根出发，做一根垂线，和曲线相交于B点，继续重复刚才的工作：
在这里插入图片描述
经过多次迭代后会越来越接近曲线的根。

上面就是牛顿迭代的几何直觉。只有将上述过程用代数表示出来，才方便我们写代码。

首先， $x_0$ 处的切线方程为： $y - y_0 = f'(x) (x - x_0)$ ，即： $f(x) = y = f'(x) (x - x_0) + f(x_0)$ 。
要求在 $x_n$ 所在切线方程交于x轴的点 $x_{n+1}$ ，即求解 $0 = f'(x) (x_{n+1} - x_n) + f(x_n)$ ，==> $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ 。
通过上面的公式可以得出： $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ ，该公式可用于求任何次方的解。
这里求的是平方根，即求解 $f(x) = x^2 = n$ 时，x的值是多少。也就是求解： $f(x) = x^2 - n = 0$ 的解。
由于 $f^{'} (x) = 2 x$ ，所以 $x_{n+1} = x_n - \frac{x_n^2 - n}{2x_n}$ ，即 $x_{n+1} = \frac{x_n + \frac{ n}{x_n}}{2}$

实现上诉公式时，一个小改进：计算机不善于做除法，所以这里修改为乘 0.5。

double sqrt_newton(double square)
{
    if (square < Min_value)
        return -1;
        
    double ret = 1.0;
    while (!equal(ret*ret, square))
    {
        ret = (ret + square / ret) * 0.5;
    }

    return ret;
}

上面的初始值是随意选择的，选择的好坏一定程度上影响了计算的速度。如，下面的方法选择的初始值是 0x5f3759df。

3、来自于Quake III源码的解法

double InvSqrt(double x)
{
    double xhalf = 0.5f * x;
    int i = *(int *)&x;             // get bits for floating value
    i = 0x5f3759df - (i >> 1);      // gives initial guess y0
    x = *(double *)&i;              // convert bits back to float
    x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
    return x;
}

4、完整代码

#include <cstdlib>
#include <cstdio>
#include <ctime>
#include <cmath>

using namespace std;

// const double Min_value = 0.0000001;
const double Min_value = 1e-7;

bool equal(double num1, double num2)
{
    if ((num1 - num2 > -Min_value) && (num1 - num2 < Min_value))
        return true;
    else
        return false;
}

double sqrt_binary_search(double square)
{
    if(square < Min_value)
        return -1;

    double low = 0.0;
    double high = square;

    // while (high - low > Min_value)
    while (!equal(high, low))
    {
        /* code */
        double mid = (high + low) / 2;
        if (mid * mid > square)
        {
            high = mid;
        }
        else
        {
            low = mid;
        }
    }

    return low;
}

double sqrt_newton(double square)
{
    if (square < Min_value)
        return -1;
        
    // double ret = 0x5f3759df; // 另一种初值
    double ret = 1.0;
    while (!equal(ret*ret, square))
    {
        ret = (ret + square / ret) * 0.5;
    }

    return ret;
}

double InvSqrt(double x)
{
    double xhalf = 0.5f * x;
    int i = *(int *)&x;             // get bits for floating value
    i = 0x5f3759df - (i >> 1);      // gives initial guess y0
    x = *(double *)&i;              // convert bits back to float
    x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
    return x;
}

int main()
{
    double x;
    scanf("%lf", &x);

    int start = clock();
    printf("sqrt_binary_search:%lf sqrt = %lf, time = %lf\n", x, 
    		sqrt_binary_search(x), (double)(clock() - start) / CLOCKS_PER_SEC);

    start = clock();
    printf("sqrt_newton:%lf sqrt = %lf, time = %lf\n", x, sqrt_newton(x), 
    		(double)(clock() - start) / CLOCKS_PER_SEC);

    start = clock();
    printf("InvSqrt:%lf sqrt = %lf, time = %lf\n", x, InvSqrt(x), 
    		(double)(clock() - start) / CLOCKS_PER_SEC);

    start = clock();
    printf("sqrt:%lf sqrt = %lf, time = %lf\n", x, sqrt(x), 
    		(double)(clock() - start) / CLOCKS_PER_SEC);
}

上述代码经过测试，可以发现 来自于Quake III源码的解法 和 cmath中的 sqrt 的速度差不多。

在这里插入图片描述

参考

旧闻一则：神秘的0x5f3759df 不可思议的Quake III源码：讲解了
FAST INVERSE SQUARE ROOT：该论文中讲解了0x5f3759df的由来。
神秘常量复出！用0x077CB531计算末尾0的个数
平方根倒数速算法 - 维基百科

Donald Su

关注

3
点赞
踩
13

收藏

觉得还不错? 一键收藏
0
评论
不用sqrt()函数，求平方根的三种方法

最近看到了这个比较有意思的题目，探究了一下。文章目录1、二分法2、牛顿法3、来自于Quake III源码的解法4、完整代码参考当然有最暴力的方法，直接遍历，(0.0, x)区间内所有的数据(也可以是 x/2)，看值是否相等，但该方法太过暴力，在此不讨论。可以采用下面的二分法、牛顿法，Quake III源码中的方法。浮点数的比较方法比较特别，不能使用平常的 x == 0的方法进行比较，下面都用到了equal比较方法，equal的代码：const double Min_value = 1e-7;.
复制链接

扫一扫