結果

問題 No.1907 DETERMINATION
ユーザー hitonanode
提出日時 2022-02-06 13:24:03
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 763 ms / 4,000 ms
コード長 5,437 bytes
コンパイル時間 1,691 ms
コンパイル使用メモリ 104,072 KB
実行使用メモリ 7,040 KB
最終ジャッジ日時 2024-12-24 12:24:32
合計ジャッジ時間 27,697 ms
ジャッジサーバーID
（参考情報） 		judge3 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 4
other AC * 63
権限があれば一括ダウンロードができます

ソースコード

diff # 

// Editorial 解
#include <cassert>
#include <iostream>
#include <utility>
#include <vector>
using namespace std;

// Library Checker Characteristic Polynomial https://judge.yosupo.jp/problem/characteristic_polynomial
// Upper Hessenberg reduction of square matrices
// Complexity: O(n^3)
// Reference:
// http://www.phys.uri.edu/nigh/NumRec/bookfpdf/f11-5.pdf
template <class Tp> void hessenberg_reduction(std::vector<std::vector<Tp>> &M) {
    assert(M.size() == M[0].size());
    const int N = M.size();
    for (int r = 0; r < N - 2; r++) {
        int piv = -1;
        for (int h = r + 1; h < N; ++h) {
            if (M[h][r] != 0) {
                piv = h;
                break;
            }
        }
        if (piv < 0) continue;
        for (int i = 0; i < N; i++) std::swap(M[r + 1][i], M[piv][i]);
        for (int i = 0; i < N; i++) std::swap(M[i][r + 1], M[i][piv]);

        const auto rinv = Tp(1) / M[r + 1][r];
        for (int i = r + 2; i < N; i++) {
            const auto n = M[i][r] * rinv;
            for (int j = 0; j < N; j++) M[i][j] -= M[r + 1][j] * n;
            for (int j = 0; j < N; j++) M[j][r + 1] += M[j][i] * n;
        }
    }
}
// Characteristic polynomial of matrix M (|xI - M|)
// Complexity: O(n^3)
// R. Rehman, I. C. Ipsen, "La Budde's Method for Computing Characteristic Polynomials," 2011.
template <class Tp> std::vector<Tp> characteristic_poly(std::vector<std::vector<Tp>> M) {
    hessenberg_reduction(M);
    const int N = M.size();
    // p[i + 1] = (Characteristic polynomial of i-th leading principal minor)
    std::vector<std::vector<Tp>> p(N + 1);
    p[0] = {1};
    for (int i = 0; i < N; i++) {
        p[i + 1].assign(i + 2, 0);
        for (int j = 0; j < i + 1; j++) p[i + 1][j + 1] += p[i][j];
        for (int j = 0; j < i + 1; j++) p[i + 1][j] -= p[i][j] * M[i][i];
        Tp betas = 1;
        for (int j = i - 1; j >= 0; j--) {
            betas *= M[j + 1][j];
            Tp hb = -M[j][i] * betas;
            for (int k = 0; k < j + 1; k++) p[i + 1][k] += hb * p[j][k];
        }
    }
    return p[N];
}
// Library Checker ここまで


template <class T>
std::vector<T> det_of_first_degree_mat(std::vector<std::vector<T>> M0, std::vector<std::vector<T>> M1) {
    const int N = M0.size();

    int multiply_by_x = 0;  // 「特定の列に x をかける」操作を行った回数
    T detAdetBinv = 1;  // 解説中の 1 / (det A det B) の値

    for (int p = 0; p < N; ++p) {
        // M1[p][p] に nonzero を持ってきて、M1 の第 p 列を掃き出す

        int pivot = -1;
        for (int row = p; row < N; ++row) {
            if (M1[row][p] != T()) {
                pivot = row;
                break;
            }
        }

        if (pivot < 0) {
            ++multiply_by_x;
            if (multiply_by_x > N) return std::vector<T>(N + 1);

            // M1 の第 p 列で pivot が見つからなかった場合、M0 + x M1 の第 p 列に x をかけたい

            // かける前に M1 の第 p 列を第 1 ~ (p - 1) 列を使って掃き出して、
            // x をかけた後で x の二次の項が出てこないようにする
            for (int row = 0; row < p; ++row) {
                T v = M1[row][p];
                M1[row][p] = 0;
                for (int i = 0; i < N; ++i) M0[i][p] -= v * M0[i][row];
            }
            for (int i = 0; i < N; ++i) swap(M0[i][p], M1[i][p]);

            --p;  // 第 p 列をもう一度やり直す この処理は高々 N 回しか走らないので全体の計算量は O(n^3) が保たれる
            continue;
        }

        if (pivot != p) {
            M1[pivot].swap(M1[p]);
            M0[pivot].swap(M0[p]);
            detAdetBinv *= -1;
        }


        // p 行目を定数倍して M1[p][p] == 1 にする
        T v = M1[p][p], vinv = v.inv();
        detAdetBinv *= v;
        for (int col = 0; col < N; ++col) {
            M0[p][col] *= vinv;
            M1[p][col] *= vinv;
        }

        // p 行目を使用して M1 の p 列目を p 行目以外ゼロにする
        for (int row = 0; row < N; ++row) {
            if (row == p) continue;
            T v = M1[row][p];
            for (int col = 0; col < N; ++col) {
                M0[row][col] -= M0[p][col] * v;
                M1[row][col] -= M1[p][col] * v;
            }
        }
    }

    // この時点で M1 = I なので det(xI + M0) を求める
    for (auto &vec : M0) {
        for (auto &x : vec) x = -x;
    }
    auto poly = characteristic_poly(M0);
    for (auto &x : poly) x *= detAdetBinv;

    poly.erase(poly.begin(), poly.begin() + multiply_by_x);
    poly.resize(N + 1);
    return poly;
}


#include <atcoder/modint>
using mint = atcoder::modint998244353;

int main() {
    cin.tie(nullptr), ios::sync_with_stdio(false);
    int N;
    cin >> N;

    vector M0(N, vector<mint>(N)), M1(N, vector<mint>(N));

    for (auto &vec : M0) {
        for (auto &x : vec) {
            int v;
            cin >> v;
            x = v;
        }
    }

    for (auto &vec : M1) {
        for (auto &x : vec) {
            int v;
            cin >> v;
            x = v;
        }
    }

    auto ret = det_of_first_degree_mat(M0, M1);
    for (auto x : ret) cout << x.val() << '\n';
}
0