Matrix/AdjugateMatrix.hpp

• View this file on GitHub
• Last update: 2023-02-11 22:28:05+08:00
• Include: #include "Matrix/AdjugateMatrix.hpp"

Depends on

• Matrix/GaussElimination.hpp
• Matrix/LinearEquation.hpp
• Template/Template.hpp

Code

#pragma once
#include "GaussElimination.hpp"
#include "LinearEquation.hpp"

template <class T>
std::vector<std::vector<T>> Adjugate_Matrix(std::vector<std::vector<T>> x) {
    int N = x.size();
    assert(N > 0);
    assert(N == (int)x[0].size());

    std::vector<std::vector<T>> m(N, std::vector<T>(2 * N));
    for (int i = 0; i < N; i++) {
        copy(begin(x[i]), end(x[i]), begin(m[i]));
        m[i][N + i] = 1;
    }
    int rank;
    auto GE = GaussElimination(m, N, true);
    rank = GE.first;
    if (rank <= N - 2) {
        return std::vector<std::vector<T>>(N, std::vector<T>(N, 0));
    }
    if (rank == N) {
        T det = GaussElimination(x).second;
        std::vector<std::vector<T>> b(N);
        for (int i = 0; i < N; i++) {
            copy(begin(m[i]) + N, end(m[i]), back_inserter(b[i]));
        }
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                b[i][j] *= det;
            }
        }
        for (int i = 0; i < N; i++) {
            for (int j = i + 1; j < N; j++) std::swap(b[i][j], b[j][i]);
        }
        return b;
    }
    auto xt = x;
    for (int i = 0; i < N; i++)
        for (int j = i + 1; j < N; j++) std::swap(xt[i][j], xt[j][i]);

    auto ps = LinearEquation<T>(xt, std::vector<T>(N, 0));
    auto qs = LinearEquation<T>(x, std::vector<T>(N, 0));
    std::vector<T> p, q;
    int r(0), c(0);
    for (int i = 0; i < ps.size(); i++) {
        for (int j = 0; j < N; j++) {
            if (ps[i][j].get() != 0) {
                p = ps[i];
                r = j;
            }
        }
    }
    for (int i = 0; i < qs.size(); i++) {
        for (int j = 0; j < N; j++) {
            if (qs[i][j].get() != 0) {
                q = qs[i];
                c = j;
            }
        }
    }

    assert(p[r].get() != 0 && q[c].get() != 0);

    std::vector<std::vector<T>> b(N, std::vector<T>(N));
    std::vector<std::vector<T>> t;
    for (int i = 0; i < N; i++) {
        if (i == r) continue;
        std::vector<T> res;
        for (int j = 0; j < N; j++)
            if (j != c) res.push_back(x[i][j]);
        t.push_back(res);
    }

    auto tdet = GaussElimination(t).second;

    b[r][c] = ((r + c) & 1) ? -tdet : tdet;
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++) {
            b[i][j] = p[i] * q[j] / (p[r] * q[c]) * b[r][c];
        }

    return b;
}

#line 2 "Matrix/GaussElimination.hpp"

#line 2 "Template/Template.hpp"

#include <bits/stdc++.h>

using i64 = std::int64_t;
#line 4 "Matrix/GaussElimination.hpp"

template <typename T>
std::pair<int, T> GaussElimination(std::vector<std::vector<T>> &a,
                                   int pivot_end = -1,
                                   bool diagonalize = false) {
    int H = a.size(), W = a[0].size();
    int rank = 0, je = pivot_end;
    if (je == -1) je = W;
    T det = 1;
    for (int j = 0; j < je; j++) {
        int idx = -1;
        for (int i = rank; i < H; i++) {
            if (a[i][j] != T(0)) {
                idx = i;
                break;
            }
        }
        if (idx == -1) {
            det = 0;
            continue;
        }
        if (rank != idx) {
            det = -det;
            std::swap(a[rank], a[idx]);
        }
        det *= a[rank][j];
        if (diagonalize && a[rank][j] != T(1)) {
            T coeff = a[rank][j].inverse();
            for (int k = j; k < W; k++) a[rank][k] *= coeff;
        }
        int is = diagonalize ? 0 : rank + 1;
        for (int i = is; i < H; i++) {
            if (i == rank) continue;
            if (a[i][j] != T(0)) {
                T coeff = a[i][j] / a[rank][j];
                for (int k = j; k < W; k++) a[i][k] -= a[rank][k] * coeff;
            }
        }
        rank++;
    }
    return make_pair(rank, det);
}
#line 3 "Matrix/LinearEquation.hpp"

template <typename T>
std::vector<std::vector<T>> LinearEquation(std::vector<std::vector<T>> a,
                                           std::vector<T> b) {
    int H = a.size(), W = a[0].size();
    for (int i = 0; i < H; i++) a[i].push_back(b[i]);
    auto p = GaussElimination(a, W, true);
    int rank = p.first;

    for (int i = rank; i < H; ++i) {
        if (a[i][W] != 0) return std::vector<std::vector<T>>{};
    }

    std::vector<std::vector<T>> res(1, std::vector<T>(W));
    std::vector<int> pivot(W, -1);
    for (int i = 0, j = 0; i < rank; ++i) {
        while (a[i][j] == 0) ++j;
        res[0][j] = a[i][W], pivot[j] = i;
    }
    for (int j = 0; j < W; ++j) {
        if (pivot[j] == -1) {
            std::vector<T> x(W);
            x[j] = 1;
            for (int k = 0; k < j; ++k) {
                if (pivot[k] != -1) x[k] = -a[pivot[k]][j];
            }
            res.push_back(x);
        }
    }
    return res;
}
#line 4 "Matrix/AdjugateMatrix.hpp"

template <class T>
std::vector<std::vector<T>> Adjugate_Matrix(std::vector<std::vector<T>> x) {
    int N = x.size();
    assert(N > 0);
    assert(N == (int)x[0].size());

    std::vector<std::vector<T>> m(N, std::vector<T>(2 * N));
    for (int i = 0; i < N; i++) {
        copy(begin(x[i]), end(x[i]), begin(m[i]));
        m[i][N + i] = 1;
    }
    int rank;
    auto GE = GaussElimination(m, N, true);
    rank = GE.first;
    if (rank <= N - 2) {
        return std::vector<std::vector<T>>(N, std::vector<T>(N, 0));
    }
    if (rank == N) {
        T det = GaussElimination(x).second;
        std::vector<std::vector<T>> b(N);
        for (int i = 0; i < N; i++) {
            copy(begin(m[i]) + N, end(m[i]), back_inserter(b[i]));
        }
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                b[i][j] *= det;
            }
        }
        for (int i = 0; i < N; i++) {
            for (int j = i + 1; j < N; j++) std::swap(b[i][j], b[j][i]);
        }
        return b;
    }
    auto xt = x;
    for (int i = 0; i < N; i++)
        for (int j = i + 1; j < N; j++) std::swap(xt[i][j], xt[j][i]);

    auto ps = LinearEquation<T>(xt, std::vector<T>(N, 0));
    auto qs = LinearEquation<T>(x, std::vector<T>(N, 0));
    std::vector<T> p, q;
    int r(0), c(0);
    for (int i = 0; i < ps.size(); i++) {
        for (int j = 0; j < N; j++) {
            if (ps[i][j].get() != 0) {
                p = ps[i];
                r = j;
            }
        }
    }
    for (int i = 0; i < qs.size(); i++) {
        for (int j = 0; j < N; j++) {
            if (qs[i][j].get() != 0) {
                q = qs[i];
                c = j;
            }
        }
    }

    assert(p[r].get() != 0 && q[c].get() != 0);

    std::vector<std::vector<T>> b(N, std::vector<T>(N));
    std::vector<std::vector<T>> t;
    for (int i = 0; i < N; i++) {
        if (i == r) continue;
        std::vector<T> res;
        for (int j = 0; j < N; j++)
            if (j != c) res.push_back(x[i][j]);
        t.push_back(res);
    }

    auto tdet = GaussElimination(t).second;

    b[r][c] = ((r + c) & 1) ? -tdet : tdet;
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++) {
            b[i][j] = p[i] * q[j] / (p[r] * q[c]) * b[r][c];
        }

    return b;
}

Back to top page