求平方根的几种方式

大家好，欢迎来到IT知识分享网。

前言

  最近在看神书《SICP》，刚看了第一章，虽然有些难啃，但感觉确实啃得确实“香”。说不上醍醐灌顶，但应该也是受益匪浅了。书中介绍了一些关于计算机数值求解的一些问题，这里抽取一点求平方根的算法，做个总结，希望可以便人便己。

一、二分法求平方根

上面说的是通用的求解方程根的方法，对于平方根来说，也是适用的。代码如下：

//sqrt with binary search  //note:l less than h double sqrtWithBiSear(double num,double l,double h){ 
     //define function f(x) = x*x -num; auto difFunc = [num](double x)->double { 
     return x*x - num; }; //make sure that there will be one root at least if(difFunc(l) *difFunc(h) > 0){ 
     return -1; } const double LIMIT = 0.00001; while(l<=h){ 
     double m = (l+h)/2; double difM = difFunc(m); if(fabs(difM) < LIMIT){ 
     return m; } double difL = difFunc(l); double difH = difFunc(h); if(difL* difM <= 0) { 
     h = m; }else if( difH * difM <= 0 ) { 
     l = m; }else { 
     return -1; } } return -1; } 

其实对于单调递增的函数来说还可以简化一点：

//sqrt with binary search  double sqrtWithBiSear1(double n){ 
     double h = max(1.0,n); double l = min(1.0,n); const double LIMIT = 0.00001; while(1){ 
     double guess = (l+h)/ 2; double guessSquare = guess*guess; double dif = guessSquare - n; if(fabs(dif) < LIMIT){ 
     return guess; } if(dif > 0){ 
     h = guess; }else { 
     l = guess; } } return -1; } 

二、牛顿法求平方根

由于很多方程其实并不存在求根公式，而且很多时候现实中的一些应用也并不需要数学意义上的精确解，为了得到方程的近似零点，牛顿等人提出了使用迭代的方式来求这个近似的零点。

double sqrtWithNewdonWay(double x, double guess){ 
     //define goodEnough? function auto goodEnough = [](double x, double guess)->bool { 
     const double LIMIT = 0.00001; if(fabs(x - guess*guess) <= LIMIT){ 
     return true; } return false; }; while(!goodEnough(x, guess)){ 
     guess = (guess + x/guess) / 2; } return guess; } 

对于牛顿法来说，它在单根附近的收敛速度较快，算法逻辑也比较简单，而且可以求代数方程的重根、复根。但是牛顿法有个问题，就是其并不是时时收敛的，和具体的函数、选取的初值都有关系。关于收敛性的讨论不在这篇博客的讨论范围之内，有兴趣的童鞋可以查阅其他资料。

三、不动点法求平方根

不动点法求解方程的根，其基本原理是将f(x) = 0转化为 一个等价的x = g(x), 或者说x(n+1) = g(x(n)) 。通过上述方程得到逐渐逼近一个解x0,使得x0 = g（x0）。则f(x0) = 0，即x0是方程f(x)的零点。

还不理解？

去看看【5】中的内容吧，里面还有说不动点的几何意义。适合小阅一番。

有一点注意的是，构造出来的迭代函数g(x)不一定收敛，需要满足一定的条件，才会收敛。具体的条件这里也就不说了（因为我也不是很懂╮(╯▽╰)╭），请各位学霸同学自行查阅资料，谢谢。

事实证明，它是收敛的。（如果构造一个g（x） = a/x ,它就不是收敛的）

实现的代码如下：

 //fixed points of function double sqrtWithFixPointMethod(double x,double guess){ 
     //define improveGuess auto improveGuess = [](double x,double guess)->double { 
     return (guess+x/guess)/2; // iterative formula }; //define goodEnough? function auto goodEnough = [](double v1, double v2)->bool { 
     const double LIMIT = 0.00001; if(fabs(v1-v2) < LIMIT) { 
     return true; } return false; }; while(!googEnough(guess,nextGuess)){ 
     guess = nextGuess; nextGuess = improveGuess(x,guess); } return nextGuess; } 

对于不动点求解方程来说，它的收敛性也是个问题。构造不同的迭代函数g(x)可能会有不同的收敛性，而且要求g(x)在区间[a,b]上的函数值也在此区间内，对给定的初值要求较高。关于不动点法收敛性更详细的讨论可以参考【6】中所述。

四、更抽象的方式

上面是三种求平方根的实现代码，在SICP中还介绍了如何用更抽象的方式来定义函数。这儿也用cpp实现了一部分，供君一览。

//abstract version of binary serach method double biSearMethod(std::function<double (double)> f,double l, double h){ 
     const double LIMIT = 0.00001; while(l<=h){ 
     double mid= (l+h)/ 2; double f_m = f(mid); if(fabs(f_m) < LIMIT){ 
     return mid; } if( f_m > 0){ 
     h = mid; }else { 
     l = mid; } } } 

 //f is transform function, guess is the first guess value double fixedPointMethod(std::function<double (double)> f, double guess){ 
     // define goodEnoughFunc with lambda expression auto goodEnoughFunc = [](double v1, double v2) -> bool { 
     const double LIMIT = 0.00001; if(fabs(v1-v2) < LIMIT) { 
     return true; } return false; }; // define improveGuessFunc with lambda expression auto improveGuessFunc = [=](double guess)->double { 
     return f(guess); }; double nextGuess = improveGuessFunc(guess); while(!goodEnoughFunc(guess,nextGuess)){ 
     guess = nextGuess; nextGuess = improveGuessFunc(guess); } return nextGuess; } 

//Derivative functions std::function<double(double)> deriva(std::function<double(double)> f) { 
     const double dx = 0.00001; return [=](double x)->double { 
     return ((f(x+dx) - f(dx)) / dx); }; } std::function<double(double)> newtonTransfrom(std::function<double(double)> g){ 
     return [=](double x)->double { 
     return x - ( g(x) / (deriva(g)) (x) ); //define:xn+1 = xn - g(x) / d(g(x)) }; } double newtonMethod(std::function<double(double)> h,double guess) { 
     return fixedPointMethod(newtonTransfrom(h),guess); } 

参考