equation.cpp

求解一元高次方程

// #pragma comment(lib,"LIBCD.lib") 
//  #include <stdlib.h>
//  #include <stdio.h>
//  #include <math.h>
//  #include <iostream.h>
//   //求方程的sturm组
//  double* qsturmz(double*P1,int*P2,int n,int* qj)
//  
//  {double  *f1,a1; 
//   double   fL=1.0;
//   int     i,j,k,j2,m;
//   f1=new double[2*(n+1)];   //存放被除式和除式的系数
//   for(i=0;i<=n;i++) *(f1+n+1+i)=*(P1+i);
//   P1=P1+n+1;                //让pl指向f‘(x)的系数
//   m=n;                
//   j=*P2;
//   P2=P2+1;
//   k=*P2;                        //被除式的次数存入j中，除式的次数存入k
//   while(fL>1e-5)
//   { 
//  	for(i=0;i<=j;i++)
//  	{
//        *(f1+i)=*(f1+n+1+i);
//        if(i<j)
//  		  *(f1+n+1+i)=*(P1+i);
//  	}                        //*被除式的系数存入*(fl+i)中，除式的系数存入*(f1+n+1+i)中    
//      *(f1+n+1+j)=0;
//      P1=P1+(m--);                  //让Pl指向余式将要存入的内存区
//      while(j>=k)
//      {
//        if(fabs(*f1)>1e-6)
//  	  {
//  		  a1=*f1;
//        for(i=0;i<=k;i++)
//        *(f1+i)=*(f1+i)-*(f1+n+1+i)*a1/(*(f1+n+1));
//  	  }
//        for(i=1;i<=j;i++)
//        *(f1+i-1)=*(f1+i);
//        *(f1+j)=0;
//         j--;
//  	}
//       while(fabs(*f1)<1e-4&&j>0)
//  	 {
//  		 for(i=1;i<=j;i++)
//           *(f1+i-1)=*(f1+i);
//           *(f1+j)=0;
//           j--;
//  	 }
//       for(fL=0,i=0;i<=j;i++)
//  	 {
//  		 *(P1+i)=-(*(f1+i));           //将余式乘一1存入pl指向的内存区
//           fL=fL+fabs(*(f1+i));
//  	 }
//  	 
//       for(i=j+1;j<=k;j++)
//       *(P1+i)=0;                  //*将P1指向内存区的乘余部分充0
//       *(++P2)=j;               //*将余式的次数j存入P2指向的内存区
//       j2=j;j=k;k=j2;
//   }
//   delete[] f1;
//   P1=P1-(++m);
//   *qj=*(--P2);            //pl指向fs(x)，qj存入fs(x)的次数
//   return P1;
//  }
//  
//  求F(x)一f(x)/fs(x)
//  void   fsturmz(double*P1,double*P4,int n,int* qj) 
//  {
//   int     i,j,k,m=n-*qj;
//   double   *f1,a1;
//   f1=new double[n+1];
//   for(i=0;i<=n;i++)
//   *(f1+i)=*(P1+i);           //将f(x)的系数存入fl指向的内存区
//   for(i=0;i<=m;i++)
//   {
//  	 a1=*f1;
//       *(P1+i)=a1/(*P4);         //*p4指向fs(x)
//       for(j=0;j<=*qj;j++)
//       *(f1+j)=*(f1+j)-(*(P4+j))*a1/(*P4);
//       for(k=1;k<=n-i;k++)
//       *(f1+k-1)=*(f1+k);
//       *(f1+n-i)=0;
//   }                       //F(x)的系数存入P1指向的内存区
//   delete[] f1;
//   for(i=m+1;i<=n;i++)
//   *(P1+i)=0;
//   for(i=0;i<=m-1;i++)    //F‘(x)的系数存入(pl+n+1)指向的内存区
//   *(P1+n+1+i)=(m-i)*(*(P1+1));
//   *(P1+n+1+m)=0;
//  }
//  
//  求多项式在x处的变号数v(x)
//  int v(double*P1,int*P2,double x,int n)         
//  {
//   int  i1,i2,i3,i4=0;
//   double *f1;
//   while(*(P2+(++i4))>0);      //求sturm组中多项式的个数
//       f1=new double[i4+1];
//   for(i1=0;i1<=i4;i1++) 
//  	 *(f1+i1)=0;
//   for(i1=0,i2=n+1;i2>=1&&*P2!=-1;i1++,P1=P1+(i2--),P2++)
//   {
//  	 for(*(f1+i1)=*P1,i3=1;i3<=*P2;i3++)   //*求sturm组中各多项式
//       *(f1+i1)=x*(*(f1+i1))+*(P1+i3);
//   }                                 //在x处的函数值*
//   for(i1=0;i1<=i4;i1++)
//   {
//       if(fabs(*(f1+i1))>1e-6)     continue;
//       for(i2=i1;i2<i4;i2++)
//       *(f1+i2)=*(f1+i2+1);
//       *(f1+i4)=0;
//       i4--;
//   }           //二去掉fo(x)，fl(x)，…，fs(x)中的0
//   for(i2=0,i1=1;i1<=i4;i1++)         //求序列fo(x)，fl(x)，…，fs(x)的变号数
//   if((*(f1+i1-1))*(*(f1+i1))<0)
//   i2=i2+1;
//   delete[] f1;
//   return  i2;
//  }
//  
//  求根的界
//  double  qgdj(double*P1,int*P2)  
//  {
//   double m;           
//   int   i;
//   m=fabs(*(P1+1));
//   for(i=1;i<*P2;i++)   
//   if(fabs(*(P1+i))>m)
//   m=fabs(*(P1+i));
//   m=1+m/fabs(*P1);
//   return m;
//  }
//  
//  求多项式在x处的函数值
//  double   dxsv(double*P1,int*P2,double x)  
//  {
//   double   f=*P1;
//   int     i;
//   for(i=1;i<=*P2;i++)
//   f=x*f+*(P1+i);
//   return  f;
//  }
//   
//  将包含根的区间分隔为几个小区间
//  void  gdg1(double*P1,int*P2,double b,int n,double* f1,int* kk,int* k2) 
//  {
//   double c;
//   int    k1=0,k3,k4;
//   
//   for(k1=1;k1<*k2;k1++)
//   {
//  	 *(kk+k1)=0;
//       *(f1+k1)=0;
//       *(f1+*k2+k1)=0;
//   }
//   *kk=*k2;
//   *f1=-b;
//   *(f1+*k2)=b;
//   while(1)
//   {
//       for(k1=0;k1<*k2;k1++)
//       if(*(kk+k1)>1)break;   //取一个根的个数大于1的区间
//       if(k1==*k2)break;
//       c=(*(f1+k1)+*(f1+*k2+k1))/2;
//       k3=v(P1,P2,c,n);
//       if(k3==v(P1,P2,*(f1+k1),n))
//           *(f1+k1)=c;
//       else if(k3==v(P1,P2,*(f1+*k2+k1),n))
//           *(f1+*k2+k1)=c;
//       else          //二将区间分隔为两个区间
//  	 {
//  		 for(k4=*k2-1;k4>1;k4--)
//           if(*(kk+k4)==0)break;       //取两个空的单元存放另一区间的端点
//           *(f1+*k2+k4)=*(f1+*k2+k1);
//           *(f1+*k2+k1)=c;
//           *(f1+k4)=c;
//           *(kk+k1)=v(P1,P2,*(f1+k1),n)-v(P1,P2,*(f1+*k2+k1),n);
//           *(kk+k4)=v(P1,P2,*(f1+k4),n)-v(P1,P2,*(f1+*k2+k4),n);
//  	 }
//   }
//  }
//  void main()
//  {
//   int    n,i,*P2,i6;
//   int    *qj=0;
//   int    *k2=0;
//   int    *kk;
//   double  *P1,fk,*P4,a,b,*f1,c;
//   double* qsturmz(double*P1,int*P2,int n,int* qj);
//   void    fsturmz(double*P1,double*P4,int n,int* qj); 
//   int     v(double*P1,int*P2,double x,int n);
//   double  qgdj(double*P1,int*P2);
//   double  dxsv(double*P1,int*P2,double x);
//   void    gdg1(double*P1,int*P2,double b,int n,double* f1,int* kk,int* k2); 
//   printf("请输入方程的次数n=");
//   scanf("%d",&n);
//   P1=new double[(n+1)*(n+2)/2]; //存入sturm组
//   P2=new int[n+1];             //存放sturm组中多项式的次数
//   for(i=0;i<=n;i++)
//   *(P2+i)=-1;
//   printf("请输入方程的%d个系数并以,号间隔\n",n+1);
//   for(i=0;i<=n;i++)
//   {
//  	 scanf("%f",&fk);
//       *(P1+i)=fk;
//   }
//   for(i=0;i<=n-1;i++)
//   *(P1+(n+1)+i)=(n-i)*(*(P1+i));  //f'(x)的系数并存入(pl+n+1)内存区
//   *P2=n;
//   *(P2+1)=n-1;                   //f(x)的次数存入*p2中
//   P4=qsturmz(P1,P2,n,qj);    //p4指向fs(x)
//   if(qj>0)                   //全局变量存入fs(x)的次数
//   {
//  	 fsturmz(P1,P4,n,qj);
//       for(i=0;i<=n;i++)
//       *(P2+i)=-1;
//       *P2=n-*qj;
//       *(P2+1)=n-*qj-1;
//       P4=qsturmz(P1,P2,n,qj);
//   }  
//   b=qgdj(P1,P2);
//   a=-b;
//   i6=v(P1,P2,a,n)-v(P1,P2,b,n); //求(a，b)内根的个数
//   if(i6==0) 
//   {
//  	 printf("此方程没有实根\n");
//  	 return;
//  	// goto loop;
//  	 exit(-1);
//   }
//   printf("此方程有%d个实根它们是:\n",i6);
//   *k2=v(P1,P2,-b,n)-v(P1,P2,b,n);     //求区间内根的个数
//   f1=new double[2*(*k2)];            //存放各小区间的端点
//   kk=new int[*k2];                   //存放各小区间内根的个数
//   gdg1(P1,P2,b,n,f1,kk,k2);         //将(a，b)分隔为i5个小区间
//   for(i=0;i<i6;i++)
//   {
//  	 c=(*(f1+i)+*(f1+i6+i))/2;
//       do
//  	 {
//  		 if(v(P1,P2,*(f1+i),n)-v(P1,P2,c,n)>0)
//  			 *(f1+i6+i)=c;
//          else  
//               *(f1+i)=c;
//          c=(*(f1+i)+*(f1+i6+i))/2;
//  	 }
//  	 while(fabs(dxsv(P1,P2,c))>1e-5);
//       printf("x%d=%9.4f\n",i+1,c);   //求各小区间内的根
//  	 delete[] f1;
//  	 delete[] kk;
//  	 delete[] P1;
//  	 delete[] P2;
//   }
//  
//  /*	loop:;*/
//  }
//  
//  {
//   int    n,i,*P2,i6;
//   int    *qj=0;
//   int    *k2=0;
//   int    *kk;
//   double  *P1,fk,*P4,a,b,*f1,c;
//   double* qsturmz(double*P1,int*P2,int n,int* qj);
//   void    fsturmz(double*P1,double*P4,int n,int* qj); 
//   int     v(double*P1,int*P2,double x,int n);
//   double  qgdj(double*P1,int*P2);
//   double  dxsv(double*P1,int*P2,double x);
//   void    gdg1(double*P1,int*P2,double b,int n,double* f1,int* kk,int* k2); 
//   printf("请输入方程的次数n=");
//   scanf("%d",&n);
//   P1=new double[(n+1)*(n+2)/2]; //存入sturm组
//   P2=new int[n+1];             //存放sturm组中多项式的次数
//   for(i=0;i<=n;i++)
//   *(P2+i)=-1;
//   printf("请输入方程的%d个系数并以,号间隔\n",n+1);
//   for(i=0;i<=n;i++)
//   {
//  	 scanf("%f",&fk);
//       *(P1+i)=fk;
//   }
//   for(i=0;i<=n-1;i++)
//   *(P1+(n+1)+i)=(n-i)*(*(P1+i));  //f'(x)的系数并存入(pl+n+1)内存区
//   *P2=n;
//   *(P2+1)=n-1;                   //f(x)的次数存入*p2中
//   P4=qsturmz(P1,P2,n,qj);    //p4指向fs(x)
//   if(qj>0)                   //全局变量存入fs(x)的次数
//   {
//  	 fsturmz(P1,P4,n,qj);
//       for(i=0;i<=n;i++)
//       *(P2+i)=-1;
//       *P2=n-*qj;
//       *(P2+1)=n-*qj-1;
//       P4=qsturmz(P1,P2,n,qj);
//   }  
//   b=qgdj(P1,P2);
//   a=-b;
//   i6=v(P1,P2,a,n)-v(P1,P2,b,n); //求(a，b)内根的个数
//   if(i6==0) 
//   {
//  	 printf("此方程没有实根\n");
//  	 return;
//  	// goto loop;
//  	 exit(-1);
//   }
//   printf("此方程有%d个实根它们是:\n",i6);
//   *k2=v(P1,P2,-b,n)-v(P1,P2,b,n);     //求区间内根的个数
//   f1=new double[2*(*k2)];            //存放各小区间的端点
//   kk=new int[*k2];                   //存放各小区间内根的个数
//   gdg1(P1,P2,b,n,f1,kk,k2);         //将(a，b)分隔为i5个小区间
//   for(i=0;i<i6;i++)
//   {
//  	 c=(*(f1+i)+*(f1+i6+i))/2;
//       do
//  	 {
//  		 if(v(P1,P2,*(f1+i),n)-v(P1,P2,c,n)>0)
//  			 *(f1+i6+i)=c;
//          else  
//               *(f1+i)=c;
//          c=(*(f1+i)+*(f1+i6+i))/2;
//  	 }
//  	 while(fabs(dxsv(P1,P2,c))>1e-5);
//       printf("x%d=%9.4f\n",i+1,c);   //求各小区间内的根
//  	 delete[] f1;
//  	 delete[] kk;
//  	 delete[] P1;
//  	 delete[] P2;
//   }
//  
//  /*	loop:;*/
//  }
//  
//  P1=new double[(n+1)*(n+2)/2]; //存入sturm组
//  P2=new int[n+1];             //存放sturm组中多项式的次数
//  for(i=0;i<=n;i++)
//  *(P2+i)=-1;
//  printf("请输入方程的%d个系数并以,号间隔\n",n+1);
//  for(i=0;i<=n;i++)
//  {
// 	 scanf("%f",&fk);
//      *(P1+i)=fk;
//  }
//  for(i=0;i<=n-1;i++)
//  *(P1+(n+1)+i)=(n-i)*(*(P1+i));  //f'(x)的系数并存入(pl+n+1)内存区
//  *P2=n;
//  *(P2+1)=n-1;                   //f(x)的次数存入*p2中
//  P4=qsturmz(P1,P2,n,qj);    //p4指向fs(x)
//  if(qj>0)                   //全局变量存入fs(x)的次数
//  {
// 	 fsturmz(P1,P4,n,qj);
//      for(i=0;i<=n;i++)
//      *(P2+i)=-1;
//      *P2=n-*qj;
//      *(P2+1)=n-*qj-1;
//      P4=qsturmz(P1,P2,n,qj);
//  }  
//  b=qgdj(P1,P2);
//  a=-b;
//  i6=v(P1,P2,a,n)-v(P1,P2,b,n); //求(a，b)内根的个数
//  if(i6==0) 
//  {
// 	 printf("此方程没有实根\n");
// 	 return;
// 	// goto loop;
// //	 exit(-1);
//  }
//  printf("此方程有%d个实根它们是:\n",i6);
//  *k2=v(P1,P2,-b,n)-v(P1,P2,b,n);     //求区间内根的个数
//  f1=new double[2*(*k2)];            //存放各小区间的端点
//  kk=new int[*k2];                   //存放各小区间内根的个数
//  gdg1(P1,P2,b,n,f1,kk,k2);         //将(a，b)分隔为i5个小区间
//  for(i=0;i<i6;i++)
//  {
// 	 c=(*(f1+i)+*(f1+i6+i))/2;
//      do
// 	 {
// 		 if(v(P1,P2,*(f1+i),n)-v(P1,P2,c,n)>0)
// 			 *(f1+i6+i)=c;
//         else  
//              *(f1+i)=c;
//         c=(*(f1+i)+*(f1+i6+i))/2;
// 	 }
// 	 while(fabs(dxsv(P1,P2,c))>1e-5);
//      printf("x%d=%9.4f\n",i+1,c);   //求各小区间内的根
// 	 delete[] f1;
// 	 delete[] kk;
// 	 delete[] P1;
// 	 delete[] P2;
//  }
// 
// /*	loop:;*/
// }