模板-线性代数
1、矩阵类
class Matrix
{
public:
int n,m;
ll a[NM][NM];
Matrix(){n=0,m=0;memset(a,0,sizeof(a));}
Matrix(int _n,int _m){n=_n,m=_m;memset(a,0,sizeof(a));}
Matrix(int _n,int _m,int** _a)
{
n=_n,m=_m;
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
a[i][j]=_a[i-1][j-1];
}
};
Matrix IdtMat(int n)
{
Matrix ans(n,n);
for(int i=1;i<=n;i++) ans.a[i][i]=1;
return ans;
}
Matrix operator+(Matrix m1,Matrix m2)
{
Matrix ans(m1.n,m1.m);
for(int i=1;i<=ans.n;i++)
for(int j=1;j<=ans.m;j++)
ans.a[i][j]=m1.a[i][j]+m2.a[i][j];
return ans;
}
Matrix operator-(Matrix m1,Matrix m2)
{
Matrix ans(m1.n,m1.m);
for(int i=1;i<=ans.n;i++)
for(int j=1;j<=ans.m;j++)
ans.a[i][j]=m1.a[i][j]-m2.a[i][j];
return ans;
}
Matrix operator*(Matrix m,ll k)
{
for(int i=1;i<=m.n;i++)
for(int j=1;j<=m.m;j++)
m.a[i][j]=(m.a[i][j]*k)%p;
return m;
}
Matrix operator*(Matrix m1,Matrix m2) //m1.m=m2.n
{
Matrix ans(m1.n,m2.m);
ll t=0;
for(int i=1;i<=m1.n;i++)
for(int k=1;k<=m1.m;k++)
if((t=m1.a[i][k]))
for(int j=1;j<=m2.m;j++)
ans.a[i][j]=(t*m2.a[k][j]%p + ans.a[i][j])%p;
return ans;
}
Matrix qpow(Matrix m,ll k)
{
Matrix ans=IdtMat(m.n);
for(;k;k>>=1){if(k&1)ans=ans*m; m=m*m;}
return ans;
}
Matrix Trn(Matrix m)
{
Matrix ans(m.n,m.m);
for(int i=1;i<=ans.n;i++)
for(int j=1;j<=ans.m;j++)
ans.a[i][j]=m.a[j][i];
return ans;
}
Matrix inv(Matrix a)
{
int n=a.n;
Matrix a2=IdtMat(n);
for(int i=1;i<=n;i++)
{
int im=0;
for(im=i;im<=n;im++)
if(a.a[im][i]) break;
if(im==n+1) return IdtMat(0);
for(int j=1;j<=n;j++)
swap(a.a[i][j],a.a[im][j]),swap(a2.a[i][j],a2.a[im][j]);
ll invi=inv(a.a[i][i]);
for(int j=1;j<=n;j++)
{
if(j==i) continue;
ll tmp=a.a[j][i]*invi%p;
for(int k=i;k<=n;k++)
a.a[j][k]=(a.a[j][k]-a.a[i][k]*tmp%p+p)%p;
for(int k=1;k<=n;k++)
a2.a[j][k]=(a2.a[j][k]-a2.a[i][k]*tmp%p+p)%p;
}
}
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a2.a[i][j]=a2.a[i][j]*inv(a.a[i][i])%p;
return a2;
}
ll det(Matrix a)
{
int n=a.n;
for(int i=1;i<=n;i++)
{
int im=0;
for(im=i;im<=n;im++)
if(a.a[im][i]) break;
if(im==n+1) return 0;
for(int j=1;j<=n;j++) swap(a.a[i][j],a.a[im][j]);
ll invi=inv(a.a[i][i]);
for(int j=1;j<=n;j++)
{
if(j==i) continue;
ll tmp=a.a[j][i]*invi%p;
for(int k=i;k<=n;k++)
a.a[j][k]=(a.a[j][k]-a.a[i][k]*tmp%p+p)%p;
}
}
ll ans=1;
for(int i=1;i<=n;i++) ans=ans*a.a[i][i]%p;
return ans;
}
2、高斯消元
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
scanf("%lf",&a[i][j]);
//枚举第i项
//1、将i~n行该项系数最大的值放在第i行
// (实际上不为零的都可以)
// 若全部系数为0,说明无唯一解
//2、由于之前的操作,第i行的前i-1项系数已清零
// 所以我们只要通过加减消元将下面几行第i项消去即可
for(int i=1;i<=n;i++)
{
int im=0;
/*for(int j=i+1;j<=n;j++)
if(fabs(a[j][i])>fabs(a[im][i])) im=j;
if(fabs(a[im][i])<1e-7){printf("No Solution\n"); return;}*/
for(im=i;im<=n;im++)
if(fabs(a[im][i])>=1e-7) break;
if(im==n+1){printf("No Solution\n"); return;}
for(int j=1;j<=n+1;j++) swap(a[i][j],a[im][j]);
for(int j=1;j<=n;j++)
{
if(j==i) continue;
double tmp=a[j][i]/a[i][i];
for(int k=i;k<=n+1;k++)
a[j][k]-=a[i][k]*tmp;
}
}
for(int i=1;i<=n;i++) a[i][n+1]/=a[i][i],a[i][i]=1;
for(int i=1;i<=n;i++)
printf("%.2f\n",a[i][n+1]);