【模板】矩阵求逆

\[PA=E, PE=P\]
按照定义，P 即为 A 的逆矩阵。矩阵乘法可以看作是对右边矩阵的一个线性变换，即：A 经过 P 的线性变换变成了 E，E 经过同样的线性变换变成了 P。因此，只需要在高斯消元 A 矩阵，将 A 变成单位矩阵的同时，维护一个单位矩阵，做与 A 完全相同的线性变换即可得到逆矩阵。

代码如下

#include <bits/stdc++.h>

using namespace std;

const int mod = 1e9 + 7;

typedef long long LL;

LL fpow(LL a, LL b) {
    LL ret = 1 % mod;
    for (; b; b >>= 1, a = a * a % mod) {
        if (b & 1) {
            ret = ret * a % mod;
        }
    }
    return ret;
}

struct matrix {
    vector<vector<LL>> mat;
    int n;
    matrix(int _n) {
        n = _n;
        mat.resize(n + 1, vector<LL>(n + 1, 0));
    }
    void identify() {
        for (int i = 1; i <= n; i++) {
            mat[i][i] = 1;
        }
    }
    vector<LL>& operator[](int x) {
        return mat[x];
    }
    void add(int i, int j, LL val) {
        for (int k = 1; k <= n; k++) {
            mat[i][k] = (mat[i][k] + mat[j][k] * val % mod + mod) % mod;
        }
    }
    void multiply(int i, LL val) {
        for (int j = 1; j <= n; j++) {
            mat[i][j] = (mat[i][j] * val % mod + mod) % mod;
        }
    }
    void print() {
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                cout << mat[i][j] << " ";
            }
            cout << endl;
        }
    }
};

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    int n;
    cin >> n;
    matrix a(n), b(n);
    b.identify();
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            cin >> a[i][j];
        }
    }
    auto gauss = [&]() -> bool {
        for (int i = 1; i <= n; i++) {
            if (a[i][i] == 0) {
                for (int j = i + 1; j <= n; j++) {
                    if (a[j][i] != 0) {
                        swap(a[i], a[j]);
                        swap(b[i], b[j]);
                        break;
                    }
                }
            }
            if (a[i][i] == 0) {
                return 0;
            }
            LL inv = fpow(a[i][i], mod - 2);
            a.multiply(i, inv);
            b.multiply(i, inv);
            for (int j = i + 1; j <= n; j++) {
                b.add(j, i, -a[j][i]);
                a.add(j, i, -a[j][i]);
            }
        }
        for (int i = n - 1; i >= 1; i--) {
            for (int j = i + 1; j <= n; j++) {
                b.add(i, j, -a[i][j]);
                a.add(i, j, -a[i][j]);
            }
        }
        return 1;
    };
    if (gauss() == 1) {
        b.print();
    } else {
        cout << "No Solution" << endl;
    }
    return 0;
}