[模板] 线性代数:矩阵/高斯消元/矩阵求逆/行列式
矩阵
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;
}
行列式
行列式的性质
行列式的显式公式
-
(展开公式)
设 \((i_1, i_2, \cdots, i_n)\) 为 \((1, 2, \cdots, n)\) 的一个排列, 共 \(n!\) 个;
\(\delta(S) = (-1)^s \in \{1, -1\}\), \(s\) 为 \(S\) 的逆序对数目.\[\left| A \right| = \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} \] -
(代数余子式)
记 \(A_{i,j}\) 表示 \(A\) 的余子式, 即去掉 \(i\) 行和 \(i\) 列剩下方阵的行列式;
\(a_{i,j}\) 表示 \(A\) 的 \(i\) 行 \(i\) 列的值.\[\left| A \right| = \sum_{i=1}^n (-1)^{i+j} a_{i,j} A_{i,j} \quad (\forall j \in \{1, 2, \cdots, n\}) \]\[\left| A \right| = \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 \right| = m_{1,1} \]即选中 \(A\) 的某一行/列, 对这一行/列的所有元素分别求余子式.
-
(消元法)
将 \(A\) 利用高斯消元将其变为上三角矩阵 \(A'\), 记 \(s\) 为高斯消元过程中进行行交换的次数, 那么
\[\left| A \right| = (-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<<'\n';
}
}
}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;
}
