#include <iostream>
#include <cmath> 
using namespace std;
double cofactor(double* detPtr,int rank,int t);   //代數(shù)余子式 
double valDet( double *detPtr, int rank);      //行列式 

template <class T>
void exchange(T& t1,T& t2){T temp;temp=t1;t1=t2;t2=temp;} //交換
 
class SquareMatrix;
class Matrix{
  public:
    friend class SquareMatrix;     //配合轉(zhuǎn)換函數(shù)食用
    Matrix(){m=n=mn=0;}        //默認構(gòu)造函數(shù) 
    Matrix(int mt,int nt);       //構(gòu)造矩陣 
    Matrix(const Matrix& mtrx);    //復制構(gòu)造函數(shù) 
    Matrix(int mt,int nt,double* a);  //數(shù)組初始化矩陣 
    
    Matrix transposeMtrx();      //轉(zhuǎn)置矩陣 
                              //初等變換
    void exchangeRow(int r1,int r2,int c=0);      //交換行 
    void multiRow(int r,int k,int c=0);        //數(shù)乘行
    void addMultiRow(int r1,int r2,int k=1,int c=0);  //r1+=k*r2
    void exchangeColumn(int c1,int c2,int r=0);      //交換列 
    void multiColumn(int c,int k,int r=0);        //數(shù)乘列 
    void addMultiColumn(int c1,int c2,int k=1,int r=0);  //c1+=k*c2

    Matrix& operator =(const Matrix& mtrx); //賦值構(gòu)造函數(shù) 
    friend istream& operator>>(istream& input,Matrix& mtrx); 
    friend ostream& operator<<(ostream& output,Matrix& mtrx); //輸出矩陣 
    friend Matrix operator*(Matrix& m1,Matrix& m2);  //矩陣乘法 
  protected:
    int m;
    int n;
    int mn;
    double* matrixPtr;
};

class SquareMatrix:public Matrix{
  public:
    SquareMatrix():Matrix(){}             //默認構(gòu)造函數(shù) 
    SquareMatrix(int mt):Matrix(mt,mt){};       //構(gòu)造函數(shù) 
    SquareMatrix(int mt,double* a):Matrix(mt,mt,a){}; //數(shù)組初始化方陣 
    SquareMatrix(const Matrix& mtrx);         //矩陣到方陣轉(zhuǎn)換 
    
    SquareMatrix transposeSqrMtrx();    //轉(zhuǎn)置方陣 
    SquareMatrix adjugateSqrMatrix();    //伴隨矩陣 
    SquareMatrix inverseSqrMatrix();    //逆矩陣 

    friend istream& operator>>(istream& input,SquareMatrix& mtrx); 
                              //輸入方陣 
    friend SquareMatrix operator *(SquareMatrix& sm1,SquareMatrix& sm2);
                              //方陣乘法 
    double getDet();                  //行列式的值 
  private:
};

Matrix::Matrix(int mt,int nt){        //初始化m*n矩陣 
  m=mt;n=nt;mn=m*n;
  matrixPtr=new double[mn];
}

Matrix::Matrix(const Matrix& mtrx){     //復制構(gòu)造函數(shù) 
  m=mtrx.m;n=mtrx.n;mn=mtrx.mn;
  matrixPtr=new double[mn];
  for(int i=0;i<mn;i++) matrixPtr[i]=mtrx.matrixPtr[i];
}

Matrix::Matrix(int mt,int nt,double* a){   //數(shù)組初始化m*n矩陣 
  m=mt;n=nt;mn=m*n;
  matrixPtr=new double[mn];
  for(int i=0;i<mn;i++)
    matrixPtr[i]=a[i];
}

istream& operator>>(istream& input,Matrix& mtrx){  //重載>> 
  if(!mtrx.m){
    cout<<"enter the m,n of matrix:";
    input>>mtrx.m>>mtrx.n;
    mtrx.mn=mtrx.m*mtrx.n;
    mtrx.matrixPtr=new double[mtrx.mn];
    cout<<"enter the matrix:"<<endl;
  }
  else cout<<"enter a "<<mtrx.m<<'*'<<mtrx.n<<" matrix:"<<endl;
  for(int i=0;i<mtrx.mn;i++) input>>mtrx.matrixPtr[i];
  return input;
}

Matrix Matrix::transposeMtrx(){           //轉(zhuǎn)置矩陣 
  Matrix mtrx(n,m);
  for(int i=0;i<n;i++)
    for(int j=0;j<m;j++)
      mtrx.matrixPtr[m*i+j]=matrixPtr[n*j+i];
  return mtrx;
}

void Matrix::exchangeRow(int r1,int r2,int c){     //交換行,默認c=0 
  for(int i=c;i<n;i++)
    exchange(matrixPtr[n*r1+i],matrixPtr[n*r2+i]);
}

void Matrix::multiRow(int r,int k,int c){       //數(shù)乘行,默認c=0 
  for(int i=c;i<n;i++)
      matrixPtr[n*r+i]*=k;
}

void Matrix::addMultiRow(int r1,int r2,int k,int c){  //r1+=k*r2,默認k=1,c=0 
  for(int i=c;i<n;i++)
    matrixPtr[n*r1+i]+=matrixPtr[n*r2+i]*k;
}

void Matrix::exchangeColumn(int c1,int c2,int r){    //交換列,默認r=0 
  for(int i=r;i<m;i++)
    exchange(matrixPtr[n*i+c1],matrixPtr[n*i+c2]);
}

void Matrix::multiColumn(int c,int k,int r){       //數(shù)乘列,默認k=1,r=0
  for(int i=r;i<m;i++)
    matrixPtr[n*i+c]*=k;
}

void Matrix::addMultiColumn(int c1,int c2,int k,int r){ //c1+=k*c2,默認r=0 
  for(int i=r;i<m;i++)
    matrixPtr[n*i+c1]+=matrixPtr[n*i+c2]*k;
}

Matrix& Matrix::operator=(const Matrix& mtrx){    //重載= 
  m=mtrx.m;n=mtrx.n;mn=m*n;
  matrixPtr=new double[mn];
  for(int i=0;i<mn;i++) matrixPtr[i]=mtrx.matrixPtr[i];
  return *this;
}

ostream& operator<<(ostream& output,Matrix& mtrx){  //重載<< 
  output<<endl;
  for(int i=0;i<mtrx.m;i++){
    for(int j=0;j<mtrx.n;j++)
      output<<mtrx.matrixPtr[mtrx.n*i+j]<<' ';
    output<<endl;
  }
  output<<endl;
  return output;
}

Matrix operator *(Matrix& m1,Matrix& m2){        //重載* 
  Matrix m3(m1.m,m2.n);
  for(int i=0;i<m3.m;i++)
    for(int j=0;j<m3.n;j++){
      double val=0;
      for(int k=0;k<m2.m;k++)
        val+=m1.matrixPtr[m1.n*i+k]*m2.matrixPtr[m2.n*k+j];
      m3.matrixPtr[m3.n*i+j]=val;
    }
  return m3;
}

//我是萌萌噠分割線-------------------------------------------------------

SquareMatrix::SquareMatrix(const Matrix& mtrx){         //構(gòu)造函數(shù) 
  m=n=mtrx.m;mn=m*n;matrixPtr=new double[mn];
  for(int i=0;i<mn;i++) matrixPtr[i]=mtrx.matrixPtr[i];
} 

istream& operator>>(istream& input,SquareMatrix& mtrx){  //重載>> 
  if(!mtrx.m){
    cout<<"enter the m of squareMatrix:";
    input>>mtrx.m;
    mtrx.n=mtrx.m;mtrx.mn=mtrx.m*mtrx.n;
    mtrx.matrixPtr=new double[mtrx.mn];
    cout<<"enter the squareMatrix:"<<endl;
  }
  else cout<<"enter a "<<mtrx.m<<" order squareMatrix:"<<endl;
  for(int i=0;i<mtrx.mn;i++) input>>mtrx.matrixPtr[i];
  return input;
}

SquareMatrix SquareMatrix::transposeSqrMtrx(){       //轉(zhuǎn)置方陣 
  return SquareMatrix((*this).transposeMtrx());
} 

SquareMatrix SquareMatrix::adjugateSqrMatrix(){      //伴隨矩陣 
  SquareMatrix aSM(m);
  for(int i=0;i<mn;i++)
    aSM.matrixPtr[i]=cofactor(matrixPtr,m,i);
  aSM=aSM.transposeSqrMtrx();
  return aSM;
}

SquareMatrix SquareMatrix::inverseSqrMatrix(){       //逆矩陣 
  double det=getDet();
  if(det==0){
    cerr<<"this is a singular matrix!"<<endl;  //判斷奇異矩陣
    return 0; 
  }
  SquareMatrix aSM(m),iSM(m);
  aSM=adjugateSqrMatrix();
  for(int i=0;i<mn;i++)
    iSM.matrixPtr[i]=aSM.matrixPtr[i]/det;
  return iSM;
}

SquareMatrix operator *(SquareMatrix& sm1,SquareMatrix& sm2){   //重載* 
  SquareMatrix sm3(sm1.m);
  for(int i=0;i<sm3.m;i++)
    for(int j=0;j<sm3.n;j++){
      double val=0;
      for(int k=0;k<sm2.m;k++)
        val+=sm1.matrixPtr[sm1.n*i+k]*sm2.matrixPtr[sm2.n*k+j];
      sm3.matrixPtr[sm3.n*i+j]=val;
    }
  return sm3;
}

double SquareMatrix::getDet(){          //行列式 
  return valDet(matrixPtr,m);
}

//又是一條萌萌噠分割線------------------------------------------

double valDet( double *detPtr, int rank)
{
  double val=0;
  if(rank==1) return detPtr[0];
  for(int i=0;i<rank;i++)          //計算余子式保存在nextDetPtr[]中 
  {
    double *nextDetPtr=new double[(rank-1)*(rank-1)];
    for(int j=0;j<rank-1;j++)
      for(int k=0;k<i;k++)
        nextDetPtr[j*(rank-1)+k]=detPtr[(j+1)*rank+k];
    for(int j=0;j<rank-1;j++)
      for(int k=i;k<rank-1;k++)
        nextDetPtr[j*(rank-1)+k]=detPtr[(j+1)*rank+k+1];
    val+=detPtr[i]*valDet(nextDetPtr,rank-1)*pow(-1.0,i);
  }
  return val;
}

double cofactor(double* detPtr,int rank,int t){           //計算代數(shù)余子式 
  double *nextDetPtr=new double[(rank-1)*(rank-1)];
  for(int i=0,j=0;i<rank*rank;i++)
    if(i>=(t/rank)*rank&&i<(t/rank)*rank+rank||!((t-i)%rank));  //如果i和t同行或同列 
    else{
      nextDetPtr[j]=detPtr[i];
      j++;
    }
  return valDet(nextDetPtr,rank-1)*pow(-1.0,t/rank+t%rank);
} 

int main(){
  cout<<endl<<"測試驅(qū)動程序-------------------"<<endl;
  /*
  cout<<endl<<"輸入任意矩陣-------------------"<<endl; 
  Matrix m1;cin>>m1;cout<<m1;
  
  cout<<endl<<"輸入任意方陣-------------------"<<endl;
  SquareMatrix sm1;cin>>sm1;cout<<sm1;
  
  cout<<endl<<"輸入3*2矩陣--------------------"<<endl;
  Matrix m2(3,2);cin>>m2;cout<<m2;
  
  cout<<endl<<"輸入2階方陣--------------------"<<endl;
  SquareMatrix sm2(2);cin>>sm2;cout<<sm2;
  */
  cout<<endl<<"數(shù)組初始化矩陣-----------------"<<endl;
  double a1[6]={1,2,3,7,8,9};
  Matrix m3(2,3,a1);cout<<m3;
  
  cout<<endl<<"數(shù)組初始化方陣-----------------"<<endl;
  double a2[4]={3,4,5,6};
  SquareMatrix sm3(2,a2);cout<<sm3; 
  
  cout<<endl<<"復制構(gòu)造方陣/矩陣--------------"<<endl;
  Matrix m4;m4=m3;Matrix m5(m3);
  cout<<m4<<m5;
  SquareMatrix sm4;sm4=sm3;SquareMatrix sm5(sm3);
  cout<<sm4<<sm5; 
  
  cout<<endl<<"矩陣/方陣乘法------------------"<<endl;
  double a3[6]={1,0,3,2,1,0},a4[9]={4,1,0,-1,1,3,2,0,1};
  Matrix m6(2,3,a3),m7(3,3,a4);
  Matrix m8=m6*m7;cout<<m8;
  
  double a5[4]={1,2,2,3},a6[4]={2,3,4,1};
  SquareMatrix sm6(2,a5),sm7(2,a6);
  SquareMatrix sm8(sm6*sm7);cout<<sm8;
  
  cout<<endl<<"矩陣轉(zhuǎn)換為方陣-----------------"<<endl;
  SquareMatrix sm9(m7);cout<<m7<<sm9;
  
  cout<<endl<<"轉(zhuǎn)置矩陣/方陣------------------"<<endl;
  Matrix m9(m6.transposeMtrx());
  cout<<m6<<m9;
  SquareMatrix sm10=sm9.transposeSqrMtrx();
  cout<<sm9<<sm10; 
  
  cout<<endl<<"初等變換-----------------------"<<endl;
  cout<<m3<<m4;
  m4.exchangeRow(0,1,2);cout<<m3<<m4;
  m4.exchangeRow(0,1);cout<<m4;
  m4.exchangeColumn(0,2);cout<<m4;
  m4.multiRow(1,2);cout<<m4;
  m4.multiColumn(1,2,1);cout<<m4;
  m4.addMultiRow(0,1);cout<<m4;
  m4.addMultiColumn(0,2,2,1);cout<<m4;
  cout<<sm3<<sm4;
  sm4.exchangeRow(0,1);cout<<sm3<<sm4;
  
  cout<<endl<<"方陣的行列式值-----------------"<<endl;
  cout<<sm3<<sm3.getDet()<<endl; 
  
  cout<<endl<<"逆矩陣-------------------------"<<endl;
  SquareMatrix sm11=sm3.inverseSqrMatrix();cout<<sm11;
  SquareMatrix sm12=sm3*sm11;cout<<sm12;
  return 0;
}