luoguP4783 [模板]矩陣求逆 線性代數
阿新 • • 發佈:2018-12-17
求\(n^2\)的矩陣的逆
翻了翻題解,看到了初等矩陣這個東西,突然想起來在看線代的時候看到過....
然後又溫習了一遍線性代數的知識
不妨設\(PA = E\),其中\(P\)是一堆初等矩陣的積(必須同時是行變換)
由於\(PA = E, PE = P\),因此\(P(A, E) = (E, P)\)
所以我們只要對矩陣\((A, E)\)來做初等變換
由於我們只做行變換
因此,兩個分塊矩陣之間互相不干擾
所以當左側的\(A\)變化為\(E\)時,右邊的\(E\)自然變成了\(P\)
複雜度\(O(n^3)\)
#include <cstdio> #include <cstring> #include <iostream> #include <algorithm> using namespace std; #define ri register int #define rep(io, st, ed) for(ri io = st; io <= ed; io ++) #define drep(io, ed, st) for(ri io = ed; io >= st; io --) #define gc getchar inline int read() { int p = 0, w = 1; char c = gc(); while(c > '9' || c < '0') { if(c == '-') w = -1; c = gc(); } while(c >= '0' && c <= '9') p = p * 10 + c - '0', c = gc(); return p * w; } const int sid = 405; const int mod = 1e9 + 7; inline void inc(int &a, int b) { a += b; if(a >= mod) a -= mod; } inline void dec(int &a, int b) { a -= b; if(a < 0) a += mod; } inline int mul(int a, int b) { return 1ll * a * b % mod; } inline int inv(int a) { int ret = 1, k = mod - 2; for( ; k; k >>= 1, a = mul(a, a)) if(k & 1) ret = mul(ret, a); return ret; } int n; int A[sid][sid], B[sid][sid]; inline int Guass() { rep(i, 1, n) { int pos = i; rep(j, i + 1, n) if(A[j][i]) pos = j; if(!A[pos][i]) return 0; swap(A[i], A[pos]); swap(B[i], B[pos]); int IA = inv(A[i][i]); rep(j, 1, n) { if(i == j) continue; int ia = mul(A[j][i], IA); rep(k, 1, n) { if(k >= i) dec(A[j][k], mul(ia, A[i][k])); dec(B[j][k], mul(ia, B[i][k])); } } } rep(i, 1, n) { int IA = inv(A[i][i]); rep(j, 1, n) B[i][j] = mul(B[i][j], IA); } return 1; } int main() { n = read(); rep(i, 1, n) rep(j, 1, n) A[i][j] = read(); rep(i, 1, n) B[i][i] = 1; if(Guass()) { rep(i, 1, n) { rep(j, 1, n) printf("%d ", B[i][j]); printf("\n"); } } else printf("No Solution\n"); return 0; }
也許下次我們可以出一道求\(AP = B\)或者\(PA = B\)的\(P\)
相信能卡死一片人QAQ