QR方法近似解一元n次方程

<table border="1" width="200" cellspacing="1" cellpadding="1"><tbody></tbody></table>#include<iostream>
#include<fstream>
#include<cmath>
#include<cstdlib>

using namespace std;

//  竟然是面向对象编程啊
class qrrt
{
      private:
              int n, Max;
              double *a, *xr, *xi, **q, eps;
      public:
             qrrt(int nn)
             {
                      int i;
                      n = nn;
                      a = new double[n+1];
                      xr = new double[n];
                      xi = new double[n];
                      q = new double*[n];
                      for(i=0; i<n; i++)
                               q[i] = new double[n];
             }

             void input();          //从文件读入多项式方程的n+1个系数
             void qr_rt();
             void qr();
             void output();
             ~qrrt()
             {
                  int i;
                  for(i=0; i<n; i++)
                  {
                         delete[] q[i];
                  }
                   delete[] q;
                   delete[] a, xr, xi;
             }
};

void qrrt::input()
{
     int i;
     char str1[20];
     ifstream fin ("xishu.txt");
     if(!fin)
     {
             cout<<"\n不能打开这个文件！"<<"xishu.txt"<<endl;

     }
     for(i=n; i>=0; i--)
              fin>>a[i];                 //读入多项式系数
     fin.close();
     Max=1000000;                            //最大迭代次数
     eps = 0.001;                   //精度
}

void qrrt::qr_rt()
{
     int i,j;
     for(i=0; i<n; i++)
              for(j=0; j<n; j++)
                       q[i][j]=0.0;
     for(j=0; j<=n-1; j++)
              q[0][j]=-a[n-j-1]/a[n];
     for(i=0; i<=n-2; i++)
              q[i+1][i]=1.0;
     qr();
}

void qrrt::qr()
{
     int m, it, i, j, k, l;
     double b,c,w,g,xy,p,qq,r,x,s,e,f,z,y;
     it=0;
     m=n;
     while(m!=0)
     {
                l=m-1;
                while((l>0) && (fabs(q[l][l-1])>eps*
                            (fabs(q[l-1][l-1])+fabs(q[l][l]))))
                            l=l-1;
                if(l==m-1)
                {
                          xr[m-1]=q[m-1][m-1];
                          xi[m-1]=0.0;
                          m=m-1;
                          it=0;
                }
                else if(l==m-2)
                {
                     b=-(q[m-1][m-1]+q[m-2][m-2]);
                     c=q[m-1][m-1]*q[m-2][m-2]-q[m-1][m-2]*q[m-2][m-1];
                     w=b*b-4.0*c;
                     y=sqrt(fabs(w));
                     if(w>0.0)
                     {
                              xy=1.0;
                              if(b<0.0)
                                       xy=-1.0;
                              xr[m-1]=(-b-xy*y)/2.0;
                              xr[m-2]=c/xr[m-1];
                              xi[m-1]=0.0;
                              xi[m-2]=0.0;
                     }
                     else
                     {
                         xr[m-1]=-b/2.0;
                         xr[m-2]=xr[m-1];
                         xi[m-1]=y/2.0;
                         xi[m-2]=-xi[m-1];
                     }
                     m=m-2;
                     it=0;
                }
                else
                {
                    if(it>=Max)
                    {
                              cout<<"\n程序停止工作"<<endl;
                              return;
                    }
                    it=it+1;
                    for(j=l+2;j<=m-1; j++)
                               q[j][j-2]=0.0;
                    for(j=l+3; j<=m-1; j++)
                               q[j][j-3]=0.0;
                    for(k=l; k<=m-2; k++)
                    {
                             if(k!=l)
                             {
                                     p=q[k][k-1];
                                     qq=q[k+1][k-1];
                                     r=0.0;
                                     if(k!=m-2)
                                               r=q[k+2][k-1];
                             }
                             else
                             {
                                 x=q[m-1][m-1]+q[m-2][m-2];
                                 y=q[m-2][m-2]*q[m-1][m-1]
                                            -q[m-2][m-1]*q[m-1][m-2];
                                 p=q[l][l]*(q[l][l]-x)+q[l][l+1]*q[l+1][l]+y;
                                 qq=q[l+1][l]*(q[l][l]+q[l+1][l+1]-x);
                                 r=q[l+1][l]*q[l+2][l];
                             }
                             if((fabs(p)+fabs(qq)+fabs(r))!=0.0)
                             {
                                  xy=1.0;
                                  if(p<0.0)
                                           xy=-1.0;
                                  s=xy*sqrt(p*p+qq*qq+r*r);
                                  if(k!=l)
                                          q[k][k-1]=-s;
                                  e=-qq/s;
                                  f=-r/s;
                                  x=-p/s;
                                  y=-x-f*r/(p+s);
                                  g=e*r/(p+s);
                                  z=-x-e*qq/(p+s);
                                  for(j=k; j<=m-1; j++)
                                  {
                                           p=x*q[k][j]+e*q[k+1][j];
                                           qq=e*q[k][j]+y*q[k+1][j];
                                           r=f*q[k][j]+g*q[k+1][j];
                                           if(k!=m-2)
                                           {
                                                     p=p+f*q[k+2][j];
                                                     qq=qq+g*q[k+2][j];
                                                     r=r+z*q[k+2][j];
                                                     q[k+2][j]=r;
                                           }
                                           q[k+1][j]=qq;
                                           q[k][j]=p;
                                  }
                                  j=k+3;
                                  if(j>=m-1)
                                            j=m-1;
                                  for(i=l; i<=j; i++)
                                  {
                                           p=x*q[i][k]+e*q[i][k+1];
                                           qq=e*q[i][k]+y*q[i][k+1];
                                           r=f*q[i][k]+g*q[i][k+1];
                                           if(k!=m-2)
                                           {
                                                     p=p+f*q[i][k+2];
                                                     qq=qq+g*q[i][k+2];
                                                     r=r+z*q[i][k+2];
                                                     q[i][k+2]=r;
                                           }
                                           q[i][k+1]=qq;
                                           q[i][k]=p;
                                  }
                             }
                    }
                }
     }
}

void qrrt::output()
{
     int i;
     char str2[20];
     cout<<"\n输出文件名:";
     cin>>str2;
     ofstream fout(str2);
     if(!fout)
     {
             cout<<"\n不能打开这个文件"<<str2<<endl;

     }
     for(i=0; i<n; i++)
     {
              fout<<xr[i]<<" + "<<xi[i]<<"i"<<endl;
              cout<<xr[i]<<" + "<<xi[i]<<"i"<<endl;
     }
     fout.close();
}

int main()
{
     qrrt root(3);
     root.input();
     root.qr_rt();
     root.output();

     return 0;
}