多项式方程实根

//2005.3 杭州

#ifndef PolyEquation_Header

#define PolyEquation_Header

 

#include <cassert>

#include <vector>

#include <cmath>

#include <algorithm>

 

#define MAXRESULT 1e10

#define MINRESULT -1e10

#define EPSILON   1e-6

 

using namespace std;

 

template<typename T>

class PolyEquation

{

public:

     PolyEquation(int _N):N(_N)

     {

         coe.reserve(N+1);

     }

     template<typename T1>

         PolyEquation(const T1 &pVal,int _N):N(_N)

     {

         coe.reserve(N+1);

         for (int i=0;i<N+1;i++)

         {

              coe.push_back(pVal[i]);

         }

     }

     T GetValue(T _Val)

     {

         T tmp1=0;

         for (int i=0;i<N+1;i++)

         {

              T tmp2=1.0;

              for (int j=0;j<i;j++)

              {

                   tmp2*=_Val;

              }

              tmp1+=coe[N-i]*tmp2;

         }

         return tmp1;

     }

     PolyEquation<T> GetDPoly()

     {

         assert(N>1);

         vector<T> tmp1;

         for (int i=0;i<N;i++)

         {

              tmp1.push_back(coe[i]*(N-i));

         }

         return PolyEquation<T>(tmp1,N-1);

     }

 

     vector<T> GetRealSolve()

     {

         vector<T> tmp1;

         if(N==0)

         {

              return tmp1;

         }

         if(N==1)

         {

 

              if(coe[0]!=0)

              {

                   tmp1.push_back(coe[1]/coe[0]);

              }

              return tmp1;

         }

         if(N==2)

         {

              if(coe[0]!=0)

              {

                   if ((coe[1]*coe[1]-coe[0]*coe[2]*4.0)>=0)

                   {

                       tmp1.push_back((-coe[1]+sqrt(coe[1]*coe[1]-coe[0]*coe[2]*4.0))/(2.0*coe[0]));

                       tmp1.push_back((-coe[1]-sqrt(coe[1]*coe[1]-coe[0]*coe[2]*4.0))/(2.0*coe[0]));

                   }

              }

              else

              {

                   vector<T> tmp2;

                   tmp2.push_back(coe[1]);

                   tmp2.push_back(coe[2]);

                   PolyEquation<T> tmp3(tmp2,1);

                   return tmp3.GetRealSolve();

              }

              return tmp1;

         }

         vector<T> tmp2=GetDPoly().GetRealSolve();

         int N2=tmp2.size();

         if (N2>0)

         {

              sort(tmp2.begin(),tmp2.end());

              DivSol(MINRESULT,tmp2[0],tmp1);

              DivSol(tmp2[N2-1],MAXRESULT,tmp1);

              if(N>1)

              {

                   for (int i=0;i<N-1;i++)

                   {

                       DivSol(tmp2[i],tmp2[i+1],tmp1);

                   }

              }

         }

         else

         {

              DivSol(MINRESULT,MAXRESULT,tmp1);

         }

         return tmp1;

     }

     void DivSol(T Val1,T Val2,vector<T> &Res)

     {

         if (fabs(Val2-Val1)<EPSILON)

         {

              Res.push_back(Val1);

              return;

         }

         T tmp1=GetValue(Val1);

         T tmp2=GetValue(Val2);

         T tmp3=GetValue((Val1+Val2)/2.0);

         if(fabs(tmp1)<EPSILON)

         {

              Res.push_back(Val1);

              return;

         }

         if(fabs(tmp2)<EPSILON)

         {

              Res.push_back(Val2);

              return;

         }

         if( (tmp1>0 && tmp2>0) || (tmp1<0 && tmp2<0) )

         {

              return;

         }

         if (fabs(tmp3)<EPSILON)

         {

              Res.push_back((Val1+Val2)/2.0);

              return;

         }

         if ( (tmp1>0 && tmp3<0) || (tmp1<0 && tmp3>0))

         {

              DivSol(Val1,(Val1+Val2)/2.0,Res);

         }

         if ( (tmp2>0 && tmp3<0) || (tmp2<0 && tmp3>0))

         {

              DivSol(Val2,(Val1+Val2)/2.0,Res);

         }

     }

private:

     int N;

     vector<T> coe;

};

 

#endif