/**
  对分法求非线性方程的根  

   算法思想:在区间[a,b]间找到所有根,从a开始搜寻,步长h.
      对每一个子区间[x_{i ,}x_i+1],    x_i+1 = x_i +h:
1.若f(xi) = 0, xi为实根,从xi +h/2 继续搜索2.若f(x_i+1) = 0, x_i+1为实根,  从  x_i+1+h/2 继续搜索
3..若f(xi) f(x_i+1) >0, 当前子区间无实根, 从x_i+1继续搜索
4...若f(xi) f(x_i+1) <0, 当前子区间有实根 ,二分搜索寻找根(当然,只能找到一个根)
注意步长h的选择, 过大导致漏根, 过小导致计算量增大.
                                                               
《数值计算方法与算法》-2 Editon -科学出版社 P86
《C#数值计算算法编程》 p207

 代码维护：2007.04.20   pengkuny
**/ 

#include<iostream>
#include<cmath>

using namespace std;

#define epsilon 0.00001  //精度

double func(double x)//方程表达式
{
    double y;
    y = (((((x-5)*x+3)*x+1)*x-7)*x+7)*x-20;  //范例方程,6次方程,至多六个实根
    return y;
}

/**//**
nNumRoots:实根个数的预估值
x[]:记录搜索到的所有实根，
xStart: 区间的左端点
xEnd :  区间的右端点
Step :  搜索求根时采用的步长
返回:求得的实根的数目
*/
int getRootBisect(int nNumRoots, double x[], double xStart, double xEnd, double Step)
{
    int n,js;
    double z,y,z1,y1,z0,y0;

    // 根的个数清0
    n = 0; 

    // 从左端点开始搜索
    z = xStart; 
    y = func(z);

    // 循环求解
    while ((z<=xEnd+Step/2.0) && (n!=nNumRoots))
    { 
        if (fabs(y)<epsilon)
        { 
            n++; 
            x[n-1] = z;
            z = z+Step/2.0; 
            y = func(z);
        }
        else
        { 
            z1 = z+Step; 
            y1 = func(z1);

            if (fabs(y1)<epsilon)
            { 
                n++; 
                x[n-1]=z1;
                z=z1+Step/2.0; 
                y=func(z);
            }
            else if (y*y1>0.0)
            { 
                y=y1; 
                z=z1;
            }
            else
            { 
                js=0;//是否找到根的标志
                while (js==0)
                { 
                    if (fabs(z1-z)<epsilon)
                    { 
                        n++; 
                        x[n-1]=(z1+z)/2.0;
                        z=z1+Step/2.0; y=func(z);
                        js=1;
                    }
                    else//对分搜索
                    { 
                        z0=(z1+z)/2.0; 
                        y0=func(z0);
                        if (fabs(y0)<epsilon)
                        { 
                            x[n]=z0; 
                            n++; 
                            js=1;
                            z=z0+Step/2.0; 
                            y=func(z);
                        }
                        else if ((y*y0)<0.0)
                        { 
                            z1=z0; 
                            y1=y0;
                        }
                        else 
                        { 
                            z=z0; 
                            y=y0;
                        }
                    }
                }
            }
        }
    }

    // 返回实根的数目
    return(n);
}

int main()
{
    cout<<"对分法求解非线性方程:"<<endl;

    double a = -2.5, b = 5;
    double *x= new double[6];

    int n = getRootBisect(6, x, a, b, 0.2);
    cout<<"共有"<<n<<"个根:"<<endl;
    for (int i =0; i<n; i++)
    {
        cout<<"  "<<x[i];
    }
    cout<<endl;

    system("pause");
    return 0;
}