求出高斯消元求解一次多项式行列式

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <iterator>
#include <limits>
#include <memory>
#include <optional>
#include <sstream>
#include <tuple>
#include <utility>
#include <vector>


template <typename Tp>
using Matrix = std::vector<std::vector<Tp>>;

template <typename Tp>
inline int width(const Matrix<Tp> &A) {
    return A.empty() ? 0 : (int)A[0].size();
}

template <typename Tp>
inline int height(const Matrix<Tp> &A) {
    return A.size();
}

template <typename Tp>
inline bool is_square_matrix(const Matrix<Tp> &A) {
    return width(A) == height(A);
}

template <typename Tp>
inline Matrix<Tp> transpose(const Matrix<Tp> &A) {
    const int w = width(A);
    const int h = height(A);
    Matrix<Tp> TA(w, std::vector<Tp>(h));
    for (int i = 0; i < h; ++i)
        for (int j = 0; j < w; ++j) TA[j][i] = A[i][j];
    return TA;
}

template <typename Tp>
inline std::vector<Tp> mat_apply(const Matrix<Tp> &A, const std::vector<Tp> &b) {
    const int w = width(A);
    const int h = height(A);
    assert((int)b.size() == w);
    std::vector<Tp> Ab(h);
    for (int i = 0; i < h; ++i)
        for (int j = 0; j < w; ++j) Ab[i] += A[i][j] * b[j];
    return Ab;
}

template <typename Tp>
inline Matrix<Tp> mat_mul(const Matrix<Tp> &A, const Matrix<Tp> &B) {
    const int wA = width(A);
    const int hA = height(A);
    assert(height(B) == wA);
    const int wB = width(B);
    Matrix<Tp> res(hA, std::vector<Tp>(wB));
    for (int i = 0; i < hA; ++i)
        for (int k = 0; k < wA; ++k)
            for (int j = 0; j < wB; ++j) res[i][j] += A[i][k] * B[k][j];
    return res;
}

template <typename Tp>
inline std::optional<Matrix<Tp>> mat_inv(Matrix<Tp> A) {
    assert(is_square_matrix(A));
    const int n = height(A);
    for (int i = 0; i < n; ++i) {
        A[i].resize(n * 2);
        A[i][n + i] = 1;
    }
    for (int i = 0; i < n; ++i) {
        int pivot = i;
        for (; pivot < n; ++pivot)
            if (A[pivot][i] != 0) break;
        if (pivot == n) return {};
        if (pivot != i) A[pivot].swap(A[i]);
        if (A[i][i] != 1) {
            const auto iv = A[i][i].inv();
            for (int j = i; j < n * 2; ++j) A[i][j] *= iv;
        }
        for (int j = i + 1; j < n; ++j) {
            const auto p = A[j][i];
            if (p == 0) continue;
            for (int k = i + 1; k < n * 2; ++k) A[j][k] -= p * A[i][k];
        }
    }
    for (int i = n - 1; i > 0; --i) {
        for (int j = i - 1; j >= 0; --j) {
            const auto p = A[j][i];
            if (p == 0) continue;
            for (int k = n; k < n * 2; ++k) A[j][k] -= p * A[i][k];
        }
    }
    for (int i = 0; i < n; ++i) A[i].erase(A[i].begin(), A[i].begin() + n);
    return A;
}

template <typename Tp>
inline Tp det(Matrix<Tp> A) {
    assert(is_square_matrix(A));
    const int n = height(A);
    Tp det      = 1;
    for (int i = 0; i < n; ++i) {
        int pivot = i;
        for (; pivot < n; ++pivot)
            if (A[pivot][i] != 0) break;
        if (pivot == n) return 0;
        if (pivot != i) {
            A[pivot].swap(A[i]);
            det = -det;
        }
        det *= A[i][i];
        const auto iv = A[i][i].inv();
        for (int j = i + 1; j < n; ++j) {
            const auto p = A[j][i] * iv;
            if (p == 0) continue;
            for (int k = i; k < n; ++k) A[j][k] -= p * A[i][k];
        }
    }
    return det;
}

template <typename Tp>
inline Matrix<Tp> to_upper_hessenberg(Matrix<Tp> A) {
    assert(is_square_matrix(A));
    const int n = height(A);
    for (int i = 0; i < n - 1; ++i) {
        int pivot = i + 1;
        for (; pivot < n; ++pivot)
            if (A[pivot][i] != 0) break;
        if (pivot == n) continue;
        if (pivot != i + 1) {
            A[pivot].swap(A[i + 1]);
            for (int j = 0; j < n; ++j) std::swap(A[j][pivot], A[j][i + 1]);
        }
        const auto iv = A[i + 1][i].inv();
        for (int j = i + 2; j < n; ++j) {
            if (A[j][i] == 0) continue;
            const auto v = A[j][i] * iv;
            for (int k = i; k < n; ++k) A[j][k] -= v * A[i + 1][k];
            for (int k = 0; k < n; ++k) A[k][i + 1] += v * A[k][j];
        }
    }
    return A;
}

// returns det(xI - A)
template <typename Tp>
inline std::vector<Tp> charpoly(const Matrix<Tp> &A) {
    const auto H = to_upper_hessenberg(A);
    const int n  = height(A);
    std::vector<std::vector<Tp>> P(n + 1);
    P[0] = {Tp(1)};
    for (int i = 1; i <= n; ++i) {
        P[i].resize(i + 1);
        for (int j = 0; j < i; ++j) P[i][j] -= H[i - 1][i - 1] * P[i - 1][j], P[i][j + 1] += P[i - 1][j];
        Tp t = 1;
        for (int j = 1; j < i; ++j) {
            t *= H[i - j][i - j - 1];
            const auto prod = t * H[i - j - 1][i - 1];
            if (prod == 0) continue;
            for (int k = 0; k < i - j; ++k) P[i][k] -= prod * P[i - j - 1][k];
        }
    }
    return P[n];
}

template <unsigned Mod>
class ModInt {
    static_assert((Mod >> 31) == 0, "`Mod` must less than 2^(31)");
    template <typename Int>
    static std::enable_if_t<std::is_integral_v<Int>, unsigned> safe_mod(Int v) {
        using D = std::common_type_t<Int, unsigned>;
        return (v %= (int)Mod) < 0 ? (D)(v + (int)Mod) : (D)v;
    }

    struct PrivateConstructor {};
    static inline PrivateConstructor private_constructor{};
    ModInt(PrivateConstructor, unsigned v) : v_(v) {}

    unsigned v_;

public:
    static unsigned mod() { return Mod; }
    static ModInt from_raw(unsigned v) { return ModInt(private_constructor, v); }
    ModInt() : v_() {}
    template <typename Int, typename std::enable_if_t<std::is_signed_v<Int>, int> = 0>
    ModInt(Int v) : v_(safe_mod(v)) {}
    template <typename Int, typename std::enable_if_t<std::is_unsigned_v<Int>, int> = 0>
    ModInt(Int v) : v_(v % Mod) {}
    unsigned val() const { return v_; }

    ModInt operator-() const { return from_raw(v_ == 0 ? v_ : Mod - v_); }
    ModInt pow(long long e) const {
        if (e < 0) return inv().pow(-e);
        for (ModInt x(*this), res(from_raw(1));; x *= x) {
            if (e & 1) res *= x;
            if ((e >>= 1) == 0) return res;
        }
    }
    ModInt inv() const {
        int x1 = 1, x3 = 0, a = val(), b = Mod;
        while (b) {
            int q = a / b, x1_old = x1, a_old = a;
            x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q;
        }
        return from_raw(x1 < 0 ? x1 + (int)Mod : x1);
    }
    template <bool Odd = (Mod & 1)>
    std::enable_if_t<Odd, ModInt> div_by_2() const {
        if (v_ & 1) return from_raw((v_ + Mod) >> 1);
        return from_raw(v_ >> 1);
    }

    ModInt &operator+=(const ModInt &a) {
        if ((v_ += a.v_) >= Mod) v_ -= Mod;
        return *this;
    }
    ModInt &operator-=(const ModInt &a) {
        if ((v_ += Mod - a.v_) >= Mod) v_ -= Mod;
        return *this;
    }
    ModInt &operator*=(const ModInt &a) {
        v_ = (unsigned long long)v_ * a.v_ % Mod;
        return *this;
    }
    ModInt &operator/=(const ModInt &a) { return *this *= a.inv(); }

    friend ModInt operator+(const ModInt &a, const ModInt &b) { return ModInt(a) += b; }
    friend ModInt operator-(const ModInt &a, const ModInt &b) { return ModInt(a) -= b; }
    friend ModInt operator*(const ModInt &a, const ModInt &b) { return ModInt(a) *= b; }
    friend ModInt operator/(const ModInt &a, const ModInt &b) { return ModInt(a) /= b; }
    friend bool operator==(const ModInt &a, const ModInt &b) { return a.v_ == b.v_; }
    friend bool operator!=(const ModInt &a, const ModInt &b) { return a.v_ != b.v_; }
    friend std::istream &operator>>(std::istream &a, ModInt &b) {
        int v;
        a >> v;
        b.v_ = safe_mod(v);
        return a;
    }
    friend std::ostream &operator<<(std::ostream &a, const ModInt &b) { return a << b.val(); }
};
#line 7 "test/matrix/adjugate_matrix.0.test.cpp"

#include <array>
#include <cassert>
#include <sstream>

template <typename Tp>
inline void Dump(const Matrix<std::array<Tp, 2>> &A) {
    const int h = height(A);
    const int w = width(A);
    for (int i = 0; i < h; ++i)
        for (int j = 0; j < w; ++j) {
            std::stringstream ss;
            ss << A[i][j][0] << " + " << A[i][j][1] << "*x";
            std::string s = ss.str();
            s.resize(30, ' ');
            std::cout << s << " \n"[j == w - 1];
        }
    std::cout << '\n';
}

// returns det(A0 + xA1)
template <typename Tp>
inline std::vector<Tp> det2(Matrix<std::array<Tp, 2>> A) {
    assert(is_square_matrix(A));

    auto sub = [](auto &a, const auto &b, Tp v, int n) {
        if (v == 0) return;
        for (int i = 0; i < n; ++i) a[i][0] -= v * b[i][0], a[i][1] -= v * b[i][1];
    };

    const int n = height(A);
    Tp m        = 1;
    for (int i = 0; i < n; ++i) {
        int pivot = i;
        for (; pivot < n; ++pivot)
            if (A[pivot][i][1] != 0) break;
        if (pivot == n) continue;
        if (pivot != i) {
            A[pivot].swap(A[i]);
            m = -m;
        }
        m *= A[i][i][1];
        const auto iv = A[i][i][1].inv();
        for (int j = 0; j < n; ++j) A[i][j][0] *= iv, A[i][j][1] *= iv;
        for (int j = 0; j < i; ++j) sub(A[j], A[i], A[j][i][1], n);
        for (int j = i + 1; j < n; ++j) sub(A[j], A[i], A[j][i][1], n);
    }
    // now set A = A0 + xA1, then A1 is an upper triangular matrix
    int t = 0;
    for (; t <= n; ++t) {
        int s = 0;
        for (; s < n; ++s)
            if (std::all_of(A[s].begin(), A[s].end(), [](const auto &a) { return a[1] == 0; })) break;
        if (s == n) break;
        for (int i = 0; i < n; ++i) A[s][i][1] = std::exchange(A[s][i][0], 0);
        for (int i = 0; i < s; ++i) sub(A[s], A[i], A[s][i][1], n);
        for (int i = s + 1; i < n; ++i) sub(A[s], A[i], A[s][i][1], n);
        if (A[s][s][1] != 0) {
            m *= A[s][s][1];
            const auto iv = A[s][s][1].inv();
            for (int j = 0; j < n; ++j) A[s][j][0] *= iv, A[s][j][1] *= iv;
            for (int i = 0; i < s; ++i) sub(A[i], A[s], A[i][s][1], n);
            for (int i = s + 1; i < n; ++i) sub(A[i], A[s], A[i][s][1], n);
        }
    }
    if (t > n) return {};
    Matrix<Tp> AA(n, std::vector<Tp>(n));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j) AA[i][j] = -A[i][j][0];
    auto res = charpoly(AA);
    res.erase(res.begin(), res.begin() + t);
    for (int i = 0; i < (int)res.size(); ++i) res[i] *= m;
    return res;
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    using mint = ModInt<998244353>;

    int n;
    std::cin >> n;

    Matrix<std::array<mint, 2>> A(n, std::vector<std::array<mint, 2>>(n));

    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j) std::cin >> A[i][j][0];
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j) std::cin >> A[i][j][1];

    const auto res = det2(A);
    for (int i = 0; i < (int)res.size(); ++i) std::cout << res[i] << ' ';

    return 0;
}