結果
問題 | No.1907 DETERMINATION |
ユーザー |
![]() |
提出日時 | 2022-02-01 00:02:20 |
言語 | C++23 (gcc 13.3.0 + boost 1.87.0) |
結果 |
AC
|
実行時間 | 1,052 ms / 4,000 ms |
コード長 | 5,292 bytes |
コンパイル時間 | 1,875 ms |
コンパイル使用メモリ | 104,496 KB |
実行使用メモリ | 5,888 KB |
最終ジャッジ日時 | 2024-12-24 12:19:06 |
合計ジャッジ時間 | 33,927 ms |
ジャッジサーバーID (参考情報) |
judge3 / judge5 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
sample | AC * 4 |
other | AC * 63 |
ソースコード
#include <cassert>#include <iostream>#include <utility>#include <vector>using namespace std;// Upper Hessenberg reduction of square matrices// Complexity: O(n^3)// Reference:// http://www.phys.uri.edu/nigh/NumRec/bookfpdf/f11-5.pdftemplate <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];}#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;}}int multiply_by_x = 0; // 基本変形の最中に M0 + M1x に x をかけた回数mint detAdetBinv = 1;for (int p = 0; p < N; ++p) {// M1[p][p] に nonzero を持ってきて、M1 の第 p 行と第 p 列を全て掃き出すint piv = -1;for (int r = p; r < N; ++r) {if (M1[r][p] != 0) {piv = r;break;}}if (piv < 0) {++multiply_by_x;if (multiply_by_x > N) break;for (int i = 0; i < N; ++i) {swap(M1[i][p], M0[i][p]);}for (int r = p - 1; r >= 0; --r) {auto v = M1[r][p];for (int i = 0; i < N; ++i) {M0[i][p] -= M0[i][r] * v;M1[i][p] -= M1[i][r] * v;}assert(M1[r][p] == 0);}--p;continue;}if (piv != p) {M1[piv].swap(M1[p]);M0[piv].swap(M0[p]);detAdetBinv *= -1;}auto v = M1[p][p], vinv = v.inv();detAdetBinv *= v;// p 行目を定数倍して M1[p][p] == 1 にするfor (int j = 0; j < N; ++j) {M0[p][j] *= vinv;M1[p][j] *= vinv;}assert(M1[p][p] == 1);// p 行目を使用して M1 の p 列目を p 行目以外ゼロにするfor (int r = 0; r < N; ++r) {if (r == p) continue;if (M1[r][p] != 0) {auto v = M1[r][p];for (int j = 0; j < N; ++j) {M0[r][j] -= M0[p][j] * v;M1[r][j] -= M1[p][j] * v;}}}// p 列目を使用して M1 の p 行目を p 列目以外ゼロにするfor (int j = p + 1; j < N; ++j) {if (M1[p][j] != 0) {auto v = M1[p][j];for (int r = 0; r < N; ++r) {M0[r][j] -= M0[r][p] * v;M1[r][j] -= M1[r][p] * v;}}}}if (multiply_by_x > N) {// 行列式がゼロであることが確定for (int i = 0; i <= N; ++i) cout << 0 << '\n';return 0;}// この時点で M1 = I なので det(x + M0) を求めるfor (auto &vec : M0) {for (auto &x : vec) x = -x;}auto poly = characteristic_poly(M0);for (auto &x : poly) x *= detAdetBinv;for (int i = 0; i < multiply_by_x; ++i) poly.erase(poly.begin());poly.resize(N + 1);for (auto a : poly) cout << a.val() << '\n';}