#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <cassert>
#include <ctime>
class MclVector
{
public:
int n;
double *Mat;
/**
type=0: 列向量
type=1: 行向量
**/
int type;
MclVector() { Mat=NULL; n=; }
MclVector(int len,double initVal=0.0)
{
n=len;
Mat=];
;i<=n;i++) Mat[i]=initVal;
type=;
}
double operator[](int id) const
{
return Mat[id];
}
double& operator[](int id)
{
return Mat[id];
}
double length() const
{
;
;i<=n;i++) sum+=Mat[i]*Mat[i];
return sqrt(sum);
}
MclVector operator*(double val) const
{
MclVector ans=MclVector(n);
;i<=n;i++) ans[i]=Mat[i]*val;
return ans;
}
MclVector operator/(double val) const
{
MclVector ans=MclVector(n);
;i<=n;i++) ans[i]=Mat[i]/val;
return ans;
}
MclVector operator+(const MclVector &newVector) const
{
MclVector ans=MclVector(n);
;i<=n;i++) ans[i]=Mat[i]+newVector[i];
return ans;
}
MclVector operator-(const MclVector &newVector) const
{
MclVector ans=MclVector(n);
;i<=n;i++) ans[i]=Mat[i]-newVector[i];
return ans;
}
MclVector operator*=(double val)
{
;i<=n;i++) Mat[i]=Mat[i]*val;
return *this;
}
MclVector operator/=(double val)
{
;i<=n;i++) Mat[i]=Mat[i]/val;
return *this;
}
MclVector operator+=(const MclVector &newVector)
{
;i<=n;i++) Mat[i]+=newVector[i];
return *this;
}
MclVector operator-=(const MclVector &newVector)
{
;i<=n;i++) Mat[i]-=newVector[i];
return *this;
}
MclVector GetTranspose() const
{
MclVector ans=*this;
ans.type=;
return ans;
}
void print() const
{
;i<=n;i++) printf("%8.3lf ",Mat[i]);
puts("");
}
};
class MclMatrix
{
public:
int row,col;
MclVector *Mat;
MclMatrix() {Mat=NULL;}
MclMatrix(int _row,int _col,double initVal=0.0)
{
row=_row;
col=_col;
Mat=];
;i<=row;i++) Mat[i]=MclVector(col,initVal);
}
void setIdentityMatrix()
{
;i<=row;i++)
{
;j<=col;j++)
{
;
;
}
}
}
MclMatrix GetTranspose() const
{
MclMatrix ans=MclMatrix(col,row);
;i<=ans.row;i++)
{
;j<=ans.col;j++)
{
ans[i][j]=Mat[j][i];
}
}
return ans;
}
void print() const
{
;i<=row;i++) Mat[i].print();
puts("");
}
MclVector& operator[](int id) const
{
return Mat[id];
}
MclVector& operator[](int id)
{
return Mat[id];
}
MclMatrix operator*(const MclMatrix &Matrix) const
{
MclMatrix ans=MclMatrix(row,Matrix.col);
;i<=row;i++)
{
;j<=Matrix.col;j++)
{
;k<=col;k++)
{
ans[i][j]+=Mat[i][k]*Matrix[k][j];
}
}
}
return ans;
}
MclMatrix operator+(const MclMatrix &Matrix) const
{
MclMatrix ans=MclMatrix(row,Matrix.col);
;i<=row;i++)
{
;j<=Matrix.col;j++)
{
ans[i][j]=Mat[i][j]+Matrix[i][j];
}
}
return ans;
}
MclMatrix operator-(const MclMatrix &Matrix) const
{
MclMatrix ans=MclMatrix(row,Matrix.col);
;i<=row;i++)
{
;j<=Matrix.col;j++)
{
ans[i][j]=Mat[i][j]-Matrix[i][j];
}
}
return ans;
}
MclVector GetCol(int colId) const
{
MclVector ans=MclVector(row);
;i<=row;i++) ans[i]=Mat[i][colId];
return ans;
}
MclVector GetRow(int rowId) const
{
MclVector ans=MclVector(row);
;i<=col;i++) ans[i]=Mat[rowId][i];
return ans;
}
MclMatrix operator*=(const MclMatrix &Matrix)
{
return *this=*this*Matrix;
}
MclMatrix operator+=(const MclMatrix &Matrix)
{
return *this=*this+Matrix;
}
MclMatrix operator-=(const MclMatrix &Matrix)
{
return *this=*this-Matrix;
}
MclMatrix operator*(double x) const
{
MclMatrix ans=*this;
;i<=row;i++)
{
;j<=col;j++)
{
ans[i][j]*=x;
}
}
return ans;
}
};
MclMatrix vectorMulVector(const MclVector &A,const MclVector& B)
{
)
{
MclMatrix ans=MclMatrix(A.n,B.n);
;i<=A.n;i++)
{
;j<=B.n;j++)
{
ans[i][j]+=A[i]*B[j];
}
}
return ans;
}
else
{
assert(A.n==B.n);
MclMatrix ans=MclMatrix(,);
;i<=A.n;i++)
{
ans[][]+=A[i]*B[i];
}
return ans;
}
}
int sgn(double x)
{
;
;
;
}
/**
将矩阵A分解为一个正交矩阵Q和一个上三角矩阵R
A为任意实数矩阵
**/
std::pair<MclMatrix,MclMatrix> QRSplit(const MclMatrix &A)
{
assert(A.col==A.row);
int n=A.row;
MclMatrix Q=MclMatrix(n,n); Q.setIdentityMatrix();
MclMatrix R=A;
;i<n;i++)
{
MclVector s=R.GetCol(i);
;j<i;j++) s[j]=;
) continue;
double c=s.length();
) c*=-sgn(R[i][i]);
MclVector u=s; u[i]-=c;
MclVector uT=s.GetTranspose();
MclMatrix H=MclMatrix(n,n);
H.setIdentityMatrix();
H=H-vectorMulVector(u,uT)*(2.0/(u.length()*u.length()));
R=H*R;
Q=Q*H;
}
return std::make_pair(Q,R);
}