/**
  對分法求非線性方程的根  

   算法思想:在區間[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;
}