[模板] 线性代数:矩阵/高斯消元/矩阵求逆/行列式

矩阵

struct tmat{
    int val[psz][psz];
    int* operator[](int p){return val[p];}
    int* operator[](int p)const{return val[p];}
    void cl(){memset(v,0,sizeof(v));}
}id{{{1,0},{0,1}}},zero{{{0,0},{0,0}}};
typedef const tmat& ctmat;

tmat operator+(ctmat a,ctmat b){
	tmat res=zero;
	rep(i,0,p-1){
		rep(j,0,p-1){
			res[i][j]=(a[i][j]+b[i][j])%nmod;
		}
	}
	return res;
}
tmat operator*(ctmat a,ctmat b){
	tmat res=zero;
	rep(i,0,p-1){
		rep(j,0,p-1){
                        if(a[i][j]==0)continue;
			rep(k,0,p-1){
				res[i][k]=(a[i][j]*b[j][k]+res[i][k])%nmod;
			}
		}
	}
	return res;
}

tmat qp(tmat a,int b){//a^b
	tmat res=id;
	while(b){
		if(b&1)res=res*a;
		a=a*a,b>>=1;
	}
	return res;
}

高斯消元

简介

求解n元一次方程组.

大概就是选 (x_i) 项系数最大的为主元, 然后用这个方程消去其他方程的 (x_i) 项的系数.

代码

// arr: n*(n+1) 矩阵
// ans: 答案
int gauss(db *ans){
    db tmp;
    rep(i,1,n){
        int p0=i;
        rep(j,i+1,n)if(fabs(a[p0][i])<fabs(a[j][i]))p0=j;
        if(fabs(a[p0][i])<eps)return 0;
        if(p0!=i)rep(j,1,n+1)swap(a[i][j],a[p0][j]);
        tmp=a[i][i];
        rep(j,i,n+1)a[i][j]/=tmp;
        rep(j,1,n){
            if(j==i)continue;
            tmp=a[j][i];
            rep(k,i,n+1)a[j][k]-=a[i][k]*tmp;
        }
    }
    rep(i,1,n)ans[i]=a[i][n+1];
    return 1;
}

行列式

行列式的性质

MIT—线性代数笔记18 行列式及其性质 - 知乎

行列式的显式公式

1. (展开公式)

   设 ((i_1, i_2, cdots, i_n)) 为 ((1, 2, cdots, n)) 的一个排列, 共 (n!) 个;
   (delta(S) = (-1)^s in {1, -1}), (s) 为 (S) 的逆序对数目.

   [left| A ight| = sum_{(i_1, i_2, cdots, i_n)}delta(i_1, i_2, cdots, i_n)a_{1,i_1}a_{2,i_2}...a_{n,i_n} ]

2. (代数余子式)

   记 (A_{i,j}) 表示 (A) 的余子式, 即去掉 (i) 行和 (i) 列剩下方阵的行列式;
   (a_{i,j}) 表示 (A) 的 (i) 行 (i) 列的值.

   [left| A ight| = sum_{i=1}^n (-1)^{i+j} a_{i,j} A_{i,j} quad (forall j in {1, 2, cdots, n}) ]
   [left| A ight| = sum_{j=1}^n (-1)^{i+j} a_{i,j} A_{i,j} quad (forall i in {1, 2, cdots, n}) ]

   递归终点: (M_1) 为 (1) 行 (1) 列的矩阵,

   [left| M_1 ight| = m_{1,1} ]

   即选中 (A) 的某一行/列, 对这一行/列的所有元素分别求余子式.

3. (消元法)

   将 (A) 利用高斯消元将其变为上三角矩阵 (A'), 记 (s) 为高斯消元过程中进行行交换的次数, 那么

   [left| A ight| = (-1)^s cdot prod_{i=1}^n A'_{i,i} ]

其中第3个方法相对易于实现, 比较常用.

代码

//利用消元法
//模意义下
//实数上类似
ll a[nsz][nsz];
ll det(int n){
	ll ans=1;
	ll v,tmp;
	rep(i,1,n){
	    int p0=i;
	    rep(j,i+1,n)if(abs(a[p0][i])<abs(a[j][i]))p0=j;
	    if(a[p0][i]==0)return 0;
	    if(p0!=i){ans=-ans;rep(j,1,n)swap(a[i][j],a[p0][j]);}
	    v=inv(a[i][i]),ans=ans*a[i][i]%nmod;
	    rep(j,1,n){
	        if(j==i)continue;
			tmp=v*a[j][i]%nmod;
	        rep(k,i,n)a[j][k]=(a[j][k]-tmp*a[i][k])%nmod;//possibly -
	    }
	}
	return getv(ans%nmod);
}

矩阵求逆

简介

对于 (n) 阶方阵 (A), 求矩阵 (A^{-1}), 满足 (A A^{-1} = ID), 其中(ID)为单位矩阵.

令矩阵 (S) 为单位矩阵. 利用高斯消元将 (A) 变为单位矩阵的同时, 对 (S) 进行相同操作.

最后的 (S = A^{-1}).

代码

struct tmat{
	ll val[nsz][nsz];
	ll *operator[](int p){return val[p];}
	const ll *operator[](int p)const{return val[p];}
	void cl(){memset(val,0,sizeof(val));}
	void getid(){
		cl();
		rep(i,1,n)val[i][i]=1;
	}
	//elementary row operations
	void swap(int x,int y){rep(i,1,n)swap(val[x][i],val[y][i]);}
	void mul(int x,ll k){rep(i,1,n)val[x][i]=getv(val[x][i]*k%nmod);}
	void add(int x,ll k,int y){rep(i,1,n)val[x][i]=getv((val[x][i]+val[y][i]*k)%nmod);}//x+=k*y
	void pr(){
		rep(i,1,n){
			rep(j,1,n)cout<<val[i][j]<<' ';
			cout<<'
';
		}
	}
}a,b;

bool invmat(tmat a,tmat &b){
	b.getid();
	rep(i,1,n){
		if(a[i][i]==0){
			rep(j,i+1,n){
				if(a[j][i]){
					a.swap(i,j),b.swap(i,j);
					break;
				}
			}
		}
		if(a[i][i]==0)return 0;
		b.mul(i,inv(a[i][i])),a.mul(i,inv(a[i][i]));
		rep(j,1,n){
			if(j==i)continue;
			b.add(j,-a[j][i],i),a.add(j,-a[j][i],i);
		}
	}
//	a.pr();
	return 1;
}