哦?今天在(luogu)上fa♂现了矩阵求逆的板子……于是就切了切。
那么我们考虑一个矩阵(A),它的逆矩阵记作(A^{-1}),其中对于矩阵这个群来讲,会有(A cdot A^{-1} = I) 其中(I)表示单位矩阵,主对角线均为(1) 。
那么我们对于矩阵(A:)
(egin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3}\ a_{2,1} & a_{2,2} & a_{2,3} \ a_{3,1} & a_{3,2} & a_{3,3} end{bmatrix} quad)
我们定义他的单位增广矩阵长底下这个样:
(egin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} & | & 1 & 0 & 0\ a_{2,1} & a_{2,2} & a_{2,3} & | & 0 & 1 & 0 \ a_{3,1} & a_{3,2} & a_{3,3} & | & 0 & 0 & 1 end{bmatrix} quad)
最终这个矩阵应该是(n imes 2m)的。还有,中间那三条线是连起来的!连起来的!但是我并不知道怎么用( m{LaTex})去搞这东西,于是(GG)
然后呢,我们对这个大矩阵进行暴力消元,直到发现左边的原矩阵变成(I)为止,右边的自然就是(A^{-1}) 那么我们思考它的正确性,以下是我想的根本不严谨只能意会的证明:
我们思考对于多个基于矩阵初等变换的变换,我们不妨称其为“变换集”。那我们思考对于原矩阵而言,从初态变成(I)的时候,经过的初等变换集虽然不唯一,但是一定组成一个变换集。那么如果右边的单位矩阵的变换集相同的话(意会证明来了),那么也就是说经历了左半个矩阵从原矩阵变成单位矩阵的全过程,那么原来的单位矩阵自然就变为了逆矩阵
好的,以上都是我不知道多长时间之前学的了(233)。啥?证明太水?没办法啊博主太弱啊……
那么思考代码:高消
(hhh)矩阵求逆·完
好吧,其实有个地方值得一说,就是我们朴素的高消都是消成上三角即可,包括求矩阵的秩,或者消元,或者求行列式。但是此处我们要求消成单位矩阵……嗯,这就很(GG).但是我从(hjc)那里听来一个很牛逼的思路:我们消到最后一定是可以把对角线消成(1),这是很简单的,但,同时如果我们再用(A_{i,i})去消同一列的其他元素,那么就一定可以将其消成(0),因为此时窝萌的(A_{i,i})是(1)啊就可以为所欲为了!那么也就是先正着消一遍,再反着消,用(1)元素去消其他的,然后就(over)了。那么大概效果图如下:
此处我们用(F**k by A_{i,i})表示被(1)元素(A_{i,i})消成(0),那么
(egin{bmatrix} 1 & F**k by A_{2,4} & F**k by A_{3,3} & F**k by A_{4,4} \ 0 & 1 & F**k by A_{3,3} &F**k by A_{4,4} \ 0 & 0& 1 & F**k by A_{4,4} \ 0 & 0 & 0 & 1end{bmatrix})
(233)此处的(F**k)是一个你不可能见过的英文单词!不可能!!
上代码吧,此处的是(LuoguP4783)
#include <cstdio>
#include <iostream>
#define MAXN 401
#define ll long long
#define mod 1000000007
using namespace std ;
ll NN, Inv ; bool flag ;
ll N, A[MAXN][MAXN << 1], i, j, k ;
inline ll qr(){
ll k = 0, f = 1 ; char c = getchar() ;
while (!isdigit(c)) {if (c == '-') f = -1 ; c = getchar() ;}
while (isdigit(c)) k = (k << 3) + (k << 1) + c - 48, c = getchar() ;
return k * f ;
}
inline ll expow(ll P, ll Q){
ll res = 1 ;
while(Q){
if (Q & 1) res = res * P % mod ;
P = P * P % mod ;
Q >>= 1 ;
} return res ;
}
inline void Gauss(){
for (i = 1 ; i <= N ; ++ i){
if (A[i][i] == 0) {flag = 1 ; return ;}
Inv = expow(A[i][i], mod - 2) , A[i][i] = 1 ;
for (j = i + 1 ; j <= NN ; ++ j) A[i][j] = A[i][j] * Inv % mod ;
for (j = 1 ; j <= N ; ++ j)
if (i == j) continue ;
else {
Inv = A[j][i], A[j][i] = 0 ;
for (k = i + 1 ; k <= NN ; ++ k){
A[j][k] = (A[j][k] - Inv * A[i][k]) % mod ;
if (A[j][k] < 0) A[j][k] += mod ;
}
}
}
}
int main(){
cin >> N ; NN = N << 1 ;
for (i = 1; i <= N ; A[i][i + N] = 1, ++ i)
for (j = 1; j <= N ; ++ j) A[i][j] = qr() ;
Gauss() ;
if (flag) {
cout << "No Solution" << endl ;
return 0 ;
}
for (i = 1; i <= N ; ++ i){
for(j = N + 1 ; j <= NN ; ++ j)
printf("%lld ", A[i][j]) ;
putchar('
') ;
}
return 0 ;
}
(Upd)出锅了出锅了哈
上面的代码忘记(swap)了,所以很(GG)……
// luogu-judger-enable-o2
#include <cstdio>
#include <iostream>
#define MAXN 401
#define ll long long
#define mod 1000000007
using namespace std ;
ll NN, Inv ; bool flag ;
ll N, A[MAXN][MAXN << 1], i, j, k ;
inline ll qr(){
ll k = 0, f = 1 ; char c = getchar() ;
while (!isdigit(c)) {if (c == '-') f = -1 ; c = getchar() ;}
while (isdigit(c)) k = (k << 3) + (k << 1) + c - 48, c = getchar() ;
return k * f ;
}
inline ll expow(ll P, ll Q){
ll res = 1 ;
while(Q){
if (Q & 1) res = res * P % mod ;
P = P * P % mod ;
Q >>= 1 ;
} return res ;
}
inline void Gauss(){
for (i = 1 ; i <= N ; ++ i){
for(j = i ; j <= N ; ++ j)
if(A[j][i]){
for(k = 1 ; k <= NN ; ++ k) swap(A[i][k], A[j][k]); break ;
}
if (A[i][i] == 0) {flag = 1 ; return ;}
Inv = expow(A[i][i], mod - 2) , A[i][i] = 1 ;
for (j = i + 1 ; j <= NN ; ++ j) A[i][j] = A[i][j] * Inv % mod ;
for (j = 1 ; j <= N ; ++ j)
if (i == j) continue ;
else {
Inv = A[j][i], A[j][i] = 0 ;
for (k = i + 1 ; k <= NN ; ++ k)
A[j][k] = ((A[j][k] - Inv * A[i][k]) % mod + mod) % mod ;
}
}
}
int main(){
cin >> N ; NN = N << 1 ;
for (i = 1; i <= N ; ++ i)
for (j = 1; j <= N ; ++ j)
A[i][j] = qr() ;
for (i = 1; i <= N ; ++ i) A[i][i + N] = 1 ;
Gauss() ;
if (flag) {
cout << "No Solution" << endl ;
return 0 ;
}
for (i = 1; i <= N ; ++ i){
for(j = N + 1 ; j <= NN ; ++ j)
printf("%lld ", A[i][j]) ;
putchar('
') ;
}
}