#737488 (C++23) No.1907 DETERMINATION

提出ソース

結果

問題	No.1907 DETERMINATION
ユーザー	👑 hitonanode
提出日時	2022-02-06 13:24:03
言語	C++23 (gcc 12.3.0 + boost 1.83.0)
結果	AC
実行時間	715 ms / 4,000 ms
コード長	5,437 bytes
コンパイル時間	1,419 ms
コンパイル使用メモリ	103,820 KB
実行使用メモリ	7,092 KB
最終ジャッジ日時	2023-08-26 00:09:39
合計ジャッジ時間	26,131 ms
ジャッジサーバーID （参考情報）	judge15 / judge14

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 4,384 KB
testcase_01	AC	2 ms 4,384 KB
testcase_02	AC	1 ms 4,380 KB
testcase_03	AC	1 ms 4,380 KB
testcase_04	AC	2 ms 4,384 KB
testcase_05	AC	2 ms 4,384 KB
testcase_06	AC	2 ms 4,380 KB
testcase_07	AC	203 ms 5,400 KB
testcase_08	AC	82 ms 4,672 KB
testcase_09	AC	139 ms 4,932 KB
testcase_10	AC	492 ms 6,744 KB
testcase_11	AC	188 ms 5,364 KB
testcase_12	AC	486 ms 6,872 KB
testcase_13	AC	458 ms 6,720 KB
testcase_14	AC	415 ms 6,840 KB
testcase_15	AC	90 ms 4,704 KB
testcase_16	AC	29 ms 4,380 KB
testcase_17	AC	462 ms 6,680 KB
testcase_18	AC	318 ms 6,000 KB
testcase_19	AC	9 ms 4,388 KB
testcase_20	AC	456 ms 6,840 KB
testcase_21	AC	42 ms 4,384 KB
testcase_22	AC	415 ms 5,488 KB
testcase_23	AC	540 ms 6,728 KB
testcase_24	AC	146 ms 5,064 KB
testcase_25	AC	1 ms 4,384 KB
testcase_26	AC	499 ms 6,856 KB
testcase_27	AC	602 ms 7,092 KB
testcase_28	AC	598 ms 6,844 KB
testcase_29	AC	593 ms 6,940 KB
testcase_30	AC	2 ms 4,380 KB
testcase_31	AC	506 ms 6,936 KB
testcase_32	AC	493 ms 6,840 KB
testcase_33	AC	494 ms 6,992 KB
testcase_34	AC	499 ms 6,928 KB
testcase_35	AC	2 ms 4,380 KB
testcase_36	AC	2 ms 4,380 KB
testcase_37	AC	2 ms 4,384 KB
testcase_38	AC	494 ms 6,940 KB
testcase_39	AC	495 ms 6,844 KB
testcase_40	AC	714 ms 6,856 KB
testcase_41	AC	493 ms 6,852 KB
testcase_42	AC	715 ms 6,848 KB
testcase_43	AC	713 ms 6,844 KB
testcase_44	AC	563 ms 6,856 KB
testcase_45	AC	603 ms 6,852 KB
testcase_46	AC	493 ms 6,964 KB
testcase_47	AC	491 ms 6,836 KB
testcase_48	AC	498 ms 6,848 KB
testcase_49	AC	525 ms 6,912 KB
testcase_50	AC	497 ms 6,904 KB
testcase_51	AC	496 ms 6,852 KB
testcase_52	AC	1 ms 4,380 KB
testcase_53	AC	508 ms 5,664 KB
testcase_54	AC	508 ms 5,752 KB
testcase_55	AC	1 ms 4,384 KB
testcase_56	AC	508 ms 5,616 KB
testcase_57	AC	508 ms 5,860 KB
testcase_58	AC	409 ms 6,844 KB
testcase_59	AC	311 ms 6,836 KB
testcase_60	AC	315 ms 6,856 KB
testcase_61	AC	400 ms 6,856 KB
testcase_62	AC	314 ms 6,848 KB
testcase_63	AC	496 ms 6,884 KB
testcase_64	AC	2 ms 4,380 KB
testcase_65	AC	2 ms 4,384 KB
testcase_66	AC	2 ms 4,380 KB

権限があれば一括ダウンロードができます

ソースコード

raw source code

// 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';
}

yukicoder

結果

テストケース

ソースコード