模板-线性代数

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]);