求给定数值的平方根

1、使用二分查找解决

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <cmath>

const double THR = 1e-4;  //结束条件

double sqr_halfsplit(double n) {
    assert(n >= 0);
    int iter = 0;
    double max = n;
    double min = 1.0;
    if (n < 1.0) {  //给定数值是否>1决定了二分的上下界
        max = 1.0;
        min = n;
    }

    double k = (min + max) / 2;
    while (1) {
        double r = k*k - n;
        printf("%d: k = %f, err = %f\n", ++iter, k, r);

        if (std::abs(r) < THR) {  //这里一定要用std::abs,即,使用cmatch中的abs;如果不加std,会使用stdlib中的abs,当给定数值<1时会得到错误结果
            return k;
        }
        else if (r > THR) {
            max = k;
        }
        else {
            min = k;
        }

        k = (min + max) / 2;
    }

    return k;
}

2、使用牛顿法解决

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <cmath>

const double THR = 1e-4;

double sqr_newton(double n) { 
    assert(n >= 0);
    double k = 1.0; 
    int iter = 0;

    double err = 0.0;
    while((err = std::abs(k*k - n)) > THR) { 
        k = (k + n/k) / 2; 
        printf("%d: k = %f, err = %f\n", ++iter, k, err);
    } 
    printf("%d: k = %f, err = %f\n", ++iter, k, err);
    
    return k; 
}

3、两种方法对比

main(int argc, char *argv[])
{
        if (argc != 2) {
                fprintf(stderr, "Usage: %s <number>\n", argv[0]);
                return -1;
        }
        double n = atof(argv[1]);

        sqr_newton(n);
        printf("---\n");
        sqr_halfsplit(n);

        return 0;
}

运行以上程序,分别传入不同的参数,可见,牛顿法的迭代次数都明显小于二分查找。而且,当传入的参数数值较大,或者THR设置较小(如1e-9)时,这种差距会更加明显

$ ./sqrt 999  
1: k = 500.000000, err = 998.000000
2: k = 250.999000, err = 249001.000000
3: k = 127.489548, err = 62001.498001
4: k = 67.662742, err = 15254.584790
5: k = 41.213573, err = 3579.246673
6: k = 32.726580, err = 699.558568
7: k = 31.626113, err = 72.029043
8: k = 31.606967, err = 1.211028
9: k = 31.606961, err = 0.000367
10: k = 31.606961, err = 0.000000
---
1: k = 500.000000, err = 249001.000000
2: k = 250.500000, err = 61751.250000
3: k = 125.750000, err = 14814.062500
4: k = 63.375000, err = 3017.390625
5: k = 32.187500, err = 37.035156
6: k = 16.593750, err = -723.647461
7: k = 24.390625, err = -404.097412
8: k = 28.289062, err = -198.728943
9: k = 30.238281, err = -84.646347
10: k = 31.212891, err = -24.755459
11: k = 31.700195, err = 5.902383
12: k = 31.456543, err = -9.485904
13: k = 31.578369, err = -1.806602
14: k = 31.639282, err = 2.044180
15: k = 31.608826, err = 0.117861
16: k = 31.593597, err = -0.844603
17: k = 31.601212, err = -0.363429
18: k = 31.605019, err = -0.122798
19: k = 31.606922, err = -0.002472
20: k = 31.607874, err = 0.057694
21: k = 31.607398, err = 0.027610
22: k = 31.607160, err = 0.012569
23: k = 31.607041, err = 0.005048
24: k = 31.606982, err = 0.001288
25: k = 31.606952, err = -0.000592
26: k = 31.606967, err = 0.000348
27: k = 31.606959, err = -0.000122
28: k = 31.606963, err = 0.000113
29: k = 31.606961, err = -0.000005

$ ./sqrt 0.625
1: k = 0.812500, err = 0.375000
2: k = 0.790865, err = 0.035156
3: k = 0.790569, err = 0.000468
4: k = 0.790569, err = 0.000000
---
1: k = 0.812500, err = 0.035156
2: k = 0.718750, err = -0.108398
3: k = 0.765625, err = -0.038818
4: k = 0.789062, err = -0.002380
5: k = 0.800781, err = 0.016251
6: k = 0.794922, err = 0.006901
7: k = 0.791992, err = 0.002252
8: k = 0.790527, err = -0.000067

4、扩展:开立方根

4.1、牛顿法的理论基础

这篇知乎的文章介绍的非常通俗易懂:如何通俗易懂地讲解牛顿迭代法求开方
这里简要概况一下原理,即为什么k = (k + n/k) / 2
求给定数值n的平方根相当于求函数y = x^2 - ny = 0时候x的解,根据牛顿法,我们随意取该函数曲线上的(k, k^2-n)做切线,斜率为2k,切线方程为\frac{y-(k^2-n)}{x-k}=2k,根据知乎的文章,我们可以把该切线看做原函数曲线的近似,同样令y = 0,则x=(k + n/k) / 2,也就是代码中的迭代公式

4.2、开立方根

套用上述过程可以解决开任意根号的问题
求给定数值n的立方根相当于求函数y = x^3 - ny = 0时候x的解
随意取该函数曲线上的(k, k^3-n)做切线,斜率为3k^2,切线方程为\frac{y-(k^3-n)}{x-k}=3k^2
y = 0,则x=(2k + n/k^2) / 3

double cuberoot_newdon(double n) {
    double k = 1.0;
    int iter = 0;

    double err = 0.0;
    while((err = std::abs(k*k*k - n)) > THR) {
        k = (k*2 + n/(k*k)) / 3;
        printf("%d: k = %f, err = %f\n", ++iter, k, err);
    }
    printf("%d: k = %f, err = %f\n", ++iter, k, err);

    return k;
}
$ ./cuberoot 1000         
1: k = 334.000000, err = 999.000000
2: k = 222.669655, err = 37258704.000000
3: k = 148.453159, err = 11039356.746705
4: k = 98.983898, err = 3270661.278411
5: k = 66.023287, err = 968825.631901
6: k = 44.091993, err = 286800.416405
7: k = 29.566121, err = 84719.414485
8: k = 20.092068, err = 24845.386955
9: k = 14.220425, err = 7110.990460
10: k = 11.128649, err = 1875.661519
11: k = 10.110596, err = 378.247997
12: k = 10.001205, err = 33.547131
13: k = 10.000000, err = 0.361651
14: k = 10.000000, err = 0.000044

$ ./cuberoot -0.7787
1: k = 0.407100, err = 1.778700
2: k = -1.294798, err = 0.846169
3: k = -1.018025, err = 1.392032
4: k = -0.929140, err = 0.276356
5: k = -0.920094, err = 0.023427
6: k = -0.920005, err = 0.000227
7: k = -0.920005, err = 0.000000
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容