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、牛顿法的理论基础
这篇知乎的文章介绍的非常通俗易懂：如何通俗易懂地讲解牛顿迭代法求开方
这里简要概况一下原理，即为什么
求给定数值n的平方根相当于求函数当
时候
的解，根据牛顿法，我们随意取该函数曲线上的
做切线，斜率为
，切线方程为
，根据知乎的文章，我们可以把该切线看做原函数曲线的近似，同样令
，则
，也就是代码中的迭代公式
4.2、开立方根
套用上述过程可以解决开任意根号的问题
求给定数值n的立方根相当于求函数当
时候
的解
随意取该函数曲线上的做切线，斜率为
，切线方程为
令，则
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
