結果

問題 No.1907 DETERMINATION
ユーザー hitonanodehitonanode
提出日時 2022-02-06 13:24:03
言語 C++23
(gcc 12.3.0 + boost 1.83.0)
結果
AC  
実行時間 693 ms / 4,000 ms
コード長 5,437 bytes
コンパイル時間 1,354 ms
コンパイル使用メモリ 104,204 KB
実行使用メモリ 7,040 KB
最終ジャッジ日時 2024-06-06 19:04:28
合計ジャッジ時間 24,656 ms
ジャッジサーバーID
（参考情報） 		judge2 / judge4
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 2 ms
5,248 KB
testcase_01 AC 1 ms
5,376 KB
testcase_02 AC 2 ms
5,376 KB
testcase_03 AC 2 ms
5,376 KB
testcase_04 AC 2 ms
5,376 KB
testcase_05 AC 2 ms
5,376 KB
testcase_06 AC 2 ms
5,376 KB
testcase_07 AC 197 ms
5,376 KB
testcase_08 AC 81 ms
5,376 KB
testcase_09 AC 135 ms
5,376 KB
testcase_10 AC 470 ms
6,912 KB
testcase_11 AC 183 ms
5,760 KB
testcase_12 AC 455 ms
6,912 KB
testcase_13 AC 440 ms
6,656 KB
testcase_14 AC 404 ms
6,784 KB
testcase_15 AC 87 ms
5,376 KB
testcase_16 AC 29 ms
5,376 KB
testcase_17 AC 444 ms
6,656 KB
testcase_18 AC 309 ms
5,760 KB
testcase_19 AC 9 ms
5,376 KB
testcase_20 AC 449 ms
6,912 KB
testcase_21 AC 41 ms
5,376 KB
testcase_22 AC 401 ms
5,760 KB
testcase_23 AC 521 ms
6,784 KB
testcase_24 AC 139 ms
5,376 KB
testcase_25 AC 2 ms
5,376 KB
testcase_26 AC 485 ms
6,912 KB
testcase_27 AC 568 ms
6,784 KB
testcase_28 AC 578 ms
6,940 KB
testcase_29 AC 570 ms
6,912 KB
testcase_30 AC 2 ms
5,376 KB
testcase_31 AC 496 ms
6,912 KB
testcase_32 AC 481 ms
6,940 KB
testcase_33 AC 487 ms
6,912 KB
testcase_34 AC 480 ms
7,040 KB
testcase_35 AC 2 ms
5,376 KB
testcase_36 AC 2 ms
5,376 KB
testcase_37 AC 2 ms
5,376 KB
testcase_38 AC 481 ms
6,784 KB
testcase_39 AC 488 ms
6,784 KB
testcase_40 AC 676 ms
6,784 KB
testcase_41 AC 481 ms
6,944 KB
testcase_42 AC 689 ms
6,912 KB
testcase_43 AC 693 ms
6,912 KB
testcase_44 AC 531 ms
6,912 KB
testcase_45 AC 594 ms
6,912 KB
testcase_46 AC 479 ms
6,912 KB
testcase_47 AC 469 ms
6,912 KB
testcase_48 AC 482 ms
6,948 KB
testcase_49 AC 497 ms
6,912 KB
testcase_50 AC 476 ms
6,912 KB
testcase_51 AC 476 ms
7,040 KB
testcase_52 AC 2 ms
5,376 KB
testcase_53 AC 495 ms
5,504 KB
testcase_54 AC 481 ms
5,760 KB
testcase_55 AC 2 ms
5,376 KB
testcase_56 AC 494 ms
5,760 KB
testcase_57 AC 498 ms
5,504 KB
testcase_58 AC 396 ms
7,040 KB
testcase_59 AC 299 ms
6,912 KB
testcase_60 AC 301 ms
6,912 KB
testcase_61 AC 387 ms
6,784 KB
testcase_62 AC 310 ms
6,912 KB
testcase_63 AC 480 ms
6,784 KB
testcase_64 AC 1 ms
5,376 KB
testcase_65 AC 2 ms
5,376 KB
testcase_66 AC 1 ms
5,376 KB
権限があれば一括ダウンロードができます

ソースコード

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