1. 程式人生 > >luoguP4783 [模板]矩陣求逆 線性代數

luoguP4783 [模板]矩陣求逆 線性代數

\(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