結果
| 問題 | No.1907 DETERMINATION |
| コンテスト | |
| ユーザー |
👑 |
| 提出日時 | 2026-01-05 23:44:16 |
| 言語 | C++23 (gcc 15.2.0 + boost 1.89.0) |
| 結果 |
AC
|
| 実行時間 | 661 ms / 4,000 ms |
| コード長 | 21,669 bytes |
| 記録 | |
| コンパイル時間 | 8,021 ms |
| コンパイル使用メモリ | 393,280 KB |
| 実行使用メモリ | 7,852 KB |
| 最終ジャッジ日時 | 2026-01-05 23:44:49 |
| 合計ジャッジ時間 | 30,723 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 63 |
ソースコード
#define USE_AC_LIBRARY
#if !defined(USE_AC_LIBRARY)
#include <bits/stdc++.h>
#elif defined(ONLINE_JUDGE)
#include <bits/stdc++.h>
#include <atcoder/all>
#else
#include <precompile.hpp>
#endif
using namespace std;
using ll = long long;
using uint = unsigned int;
using ull = unsigned long long;
namespace atcoder {};
using namespace atcoder;
namespace suisen {};
using namespace suisen;
// cp-library-cpp
#if defined(ATCODER_MODINT_HPP)
using mint = modint998244353;
#endif
#define rep_kind(a, b, c, d, e, ...) e
#define rep1(i, r) for (int i = 0; i < (int)(r); ++i)
#define rep2(i, l, r) for (int i = (l); i < (int)(r); ++i)
#define rep3(i, l, r, d) for (int i = (l); i < (int)(r); i += (d))
#define rep_r1(i, r) for (int i = (r) - 1; i >= 0; --i)
#define rep_r2(i, l, r) for (int i = (r) - 1; i >= (l); --i)
#define rep_r3(i, l, r, d) for (int i = (r) - 1; i >= (l); i -= (d))
#define rep(...) rep_kind(__VA_ARGS__, rep3, rep2, rep1)(__VA_ARGS__)
#define rep_r(...) rep_kind(__VA_ARGS__, rep_r3, rep_r2, rep_r1)(__VA_ARGS__)
#define rep_t_kind(a, b, c, d, e, f, ...) f
#define rep_t1(type, i, r) for (type i = 0; i < (type)(r); ++i)
#define rep_t2(type, i, l, r) for (type i = (l); i < (type)(r); ++i)
#define rep_t3(type, i, l, r, d) for (type i = (l); i < (type)(r); i += (d))
#define rep_t_r1(type, i, r) for (type i = (r) - 1; i >= 0; --i)
#define rep_t_r2(type, i, l, r) for (type i = (r) - 1; i >= (l); --i)
#define rep_t_r3(type, i, l, r, d) for (type i = (r) - 1; i >= (l); i -= (d))
#define rep_t(...) rep_t_kind(__VA_ARGS__, rep_t3, rep_t2, rep_t1)(__VA_ARGS__)
#define rep_t_r(...) \
rep_t_kind(__VA_ARGS__, rep_t_r3, rep_t_r2, rep_t_r1)(__VA_ARGS__)
#if defined(ATCODER_MODINT_HPP)
template <int M> istream &operator>>(istream &is, static_modint<M> &x) {
ll val;
is >> val;
x = val;
return is;
}
template <int M> ostream &operator<<(ostream &os, const static_modint<M> &x) {
os << x.val();
return os;
}
#endif
inline int slot_index() {
static int idx = ios_base::xalloc();
return idx;
}
struct separator {
string s;
};
inline separator sep(string s) { return {std::move(s)}; }
inline void cleanup(ios_base::event ev, ios_base &ios, int idx) {
if (ev == ios_base::erase_event || ev == ios_base::copyfmt_event) {
void *&slot = ios.pword(idx);
delete static_cast<string *>(slot);
slot = nullptr;
}
}
inline ostream &operator<<(ostream &os, const separator &m) {
void *&slot = os.pword(slot_index());
if (!slot) {
slot = new string(m.s);
os.register_callback(&cleanup, slot_index());
} else {
*static_cast<string *>(slot) = m.s;
}
return os;
}
inline const char *current_sep(ios_base &ios) {
if (auto p = static_cast<string *>(ios.pword(slot_index())); p) {
return p->c_str();
}
return " ";
}
template <class T1, class T2>
istream &operator>>(istream &is, pair<T1, T2> &p) {
is >> p.first >> p.second;
return is;
}
template <class T1, class T2>
ostream &operator<<(ostream &os, const pair<T1, T2> &p) {
os << p.first << current_sep(os) << p.second;
return os;
}
template <class T1, class T2, class T3>
istream &operator>>(istream &is, tuple<T1, T2, T3> &t3) {
is >> get<0>(t3) >> get<1>(t3) >> get<2>(t3);
return is;
}
template <class T1, class T2, class T3>
ostream &operator<<(ostream &os, const tuple<T1, T2, T3> &t3) {
os << get<0>(t3) << current_sep(os) << get<1>(t3) << current_sep(os)
<< get<2>(t3);
return os;
}
template <class T>
concept char_range = ranges::input_range<T> &&
same_as<remove_cv_t<ranges::range_value_t<T>>, char>;
template <class T>
concept nonstring_range =
(ranges::input_range<T> || ranges::input_range<const T>) &&
!is_convertible_v<T, string_view> && !char_range<T> && !char_range<const T>;
template <nonstring_range T> istream &operator>>(istream &is, T &r) {
for (auto &e : r) {
is >> e;
}
return is;
}
template <nonstring_range T> ostream &operator<<(ostream &os, const T &r) {
bool is_first = true;
for (auto &e : r) {
if (is_first) {
is_first = false;
} else {
os << current_sep(os);
}
os << e;
}
return os;
}
template <class... T> void in(T &...args) { (cin >> ... >> args); }
template <class... T> void in_z(T &...args) {
(cin >> ... >> args);
(..., --args);
}
#define in_d(type, ...) \
type __VA_ARGS__; \
in(__VA_ARGS__)
#define in_dz(type, ...) \
type __VA_ARGS__; \
in_z(__VA_ARGS__)
template <bool do_flush = false, class Hd, class... Tl>
void out(const Hd &hd, const Tl &...tl) {
cout << hd;
if constexpr (sizeof...(Tl)) {
(cout << ... << (cout << current_sep(cout), tl));
}
cout << endl;
if (do_flush) {
cout << flush;
}
}
template <bool do_flush = false, class Hd, class... Tl>
void err(const Hd &hd, const Tl &...tl) {
cerr << hd;
if constexpr (sizeof...(Tl)) {
(cerr << ... << (cerr << current_sep(cerr), tl));
}
cerr << endl;
if (do_flush) {
cerr << flush;
}
}
auto &out_sep(string s) { return cout << sep(s); }
auto &err_sep(string s) { return cerr << sep(s); }
#define dir(dx, dy) for (auto [dx, dy] : {pair{1, 0}, {0, 1}, {-1, 0}, {0, -1}})
#define dir_8(dx, dy) \
for (auto [dx, dy] : {pair{1, 0}, \
{1, 1}, \
{0, 1}, \
{-1, 1}, \
{-1, 0}, \
{-1, -1}, \
{0, -1}, \
{1, -1}})
#define all(v) (v).begin(), (v).end()
#define all_r(v) (v).rbegin(), (v).rend()
#define iter(v, l, r) ((v).begin() + l), ((v).begin() + r)
#define dedup(v) (v).erase(unique(all(v)), (v).end())
#define yn(b) ((b) ? "Yes" : "No")
#define yes() yn(true)
#define no() yn(false)
template <class T1, class T2> bool chmin(T1 &l, const T2 &r) {
if (r < l) {
l = r;
return true;
}
return false;
}
template <class T1, class T2> bool chmax(T1 &l, const T2 &r) {
if (r > l) {
l = r;
return true;
}
return false;
}
template <class T> void increment(T &v) {
for (auto &e : v) {
++e;
}
}
template <class T> void decrement(T &v) {
for (auto &e : v) {
--e;
}
}
template <class It> iterator_traits<It>::value_type sum_of(It first, It last) {
using T = iterator_traits<It>::value_type;
return accumulate(first, last, T{});
}
template <class It, class U>
iterator_traits<It>::value_type sum_of(It first, It last, U init) {
using T = iterator_traits<It>::value_type;
return accumulate(first, last, T(init));
}
template <class Hd1, class Hd2, class... Tl>
auto min_of(Hd1 hd1, Hd2 hd2, Tl... tl) {
if constexpr (sizeof...(Tl) == 0) {
return min(hd1, hd2);
} else {
return min(hd1, min_of(hd2, tl...));
}
}
template <class Hd1, class Hd2, class... Tl>
auto max_of(Hd1 hd1, Hd2 hd2, Tl... tl) {
if constexpr (sizeof...(Tl) == 0) {
return max(hd1, hd2);
} else {
return max(hd1, max_of(hd2, tl...));
}
}
template <class Hd1, class Hd2, class... Tl>
auto gcd_of(Hd1 hd1, Hd2 hd2, Tl... tl) {
if constexpr (sizeof...(Tl) == 0) {
return gcd(hd1, hd2);
} else {
return gcd(hd1, gcd_of(hd2, tl...));
}
}
template <class Hd1, class Hd2, class... Tl>
auto lcm_of(Hd1 hd1, Hd2 hd2, Tl... tl) {
if constexpr (sizeof...(Tl) == 0) {
return lcm(hd1, hd2);
} else {
return lcm(hd1, lcm_of(hd2, tl...));
}
}
template <class T, size_t N>
auto make_vector(vector<size_t> &sizes, T const &x) {
if constexpr (N == 1) {
return vector(sizes[0], x);
} else {
size_t size = sizes[N - 1];
sizes.pop_back();
return vector(size, make_vector<T, N - 1>(sizes, x));
}
}
template <class T, size_t N>
auto make_vector(size_t const (&sizes)[N], T const &x = T()) {
vector<size_t> s(N);
rep(i, N) { s[i] = sizes[N - i - 1]; }
return make_vector<T, N>(s, x);
}
template <class T, size_t N> auto make_vector(vector<int> &sizes, T const &x) {
if constexpr (N == 1) {
return vector(sizes[0], x);
} else {
int size = sizes[N - 1];
sizes.pop_back();
return vector(size, make_vector<T, N - 1>(sizes, x));
}
}
template <class T, size_t N>
auto make_vector(int const (&sizes)[N], T const &x = T()) {
vector<int> s(N);
rep(i, N) { s[i] = sizes[N - i - 1]; }
return make_vector<T, N>(s, x);
}
template <class T, size_t Hd, size_t... Tl> auto make_array() {
if constexpr (sizeof...(Tl) == 0) {
return array<T, Hd>{};
} else {
return array<decltype(make_array<T, Tl...>()), Hd>{};
}
}
constexpr int inf32 = 1'000'000'001;
constexpr ll inf64 = 2'000'000'000'000'000'001;
constexpr ll ten_powers[19] = {1,
10,
100,
1000,
10000,
100000,
1000000,
10000000,
100000000,
1000000000,
10000000000,
100000000000,
1000000000000,
10000000000000,
100000000000000,
1000000000000000,
10000000000000000,
100000000000000000,
1000000000000000000};
constexpr ll two_powers[63] = {1,
2,
4,
8,
16,
32,
64,
128,
256,
512,
1024,
2048,
4096,
8192,
16384,
32768,
65536,
131072,
262144,
524288,
1048576,
2097152,
4194304,
8388608,
16777216,
33554432,
67108864,
134217728,
268435456,
536870912,
1073741824,
2147483648,
4294967296,
8589934592,
17179869184,
34359738368,
68719476736,
137438953472,
274877906944,
549755813888,
1099511627776,
2199023255552,
4398046511104,
8796093022208,
17592186044416,
35184372088832,
70368744177664,
140737488355328,
281474976710656,
562949953421312,
1125899906842624,
2251799813685248,
4503599627370496,
9007199254740992,
18014398509481984,
36028797018963968,
72057594037927936,
144115188075855872,
288230376151711744,
576460752303423488,
1152921504606846976,
2305843009213693952,
4611686018427387904};
template <class T> using pq = priority_queue<T, vector<T>, greater<>>;
template <class S, S e> constexpr S e_const() { return e; }
template <class S1, S1 e1, class S2, S2 e2> constexpr pair<S1, S2> e_pair() {
return {e1, e2};
}
template <class S> S op_min(S a, S b) { return min(a, b); }
template <class S> S op_max(S a, S b) { return max(a, b); }
template <class S> S op_add(S a, S b) { return a + b; }
template <class S1, class S2>
pair<S1, S2> op_add_pair(pair<S1, S2> a, pair<S1, S2> b) {
auto [a1, a2] = a;
auto [b1, b2] = b;
return {a1 + b1, a2 + b2};
}
template <class F1, class F2, class S> S mapping_affine(pair<F1, F2> f, S x) {
auto [a, b] = f;
return a * x + b;
}
template <class F1, class F2, class S1, class S2>
pair<S1, S2> mapping_len_and_affine(pair<F1, F2> f, pair<S1, S2> x) {
auto [a, b] = f;
auto [x1, x2] = x;
return {x1, a * x2 + b * x1};
}
template <class F1, class F2>
pair<F1, F2> composition_affine(pair<F1, F2> f, pair<F1, F2> g) {
auto [a_f, b_f] = f;
auto [a_g, b_g] = g;
return {a_f * a_g, a_f * b_g + b_f};
}
template <class M, class S, S e> constexpr M mint_const() { return M(e); }
template <class M1, class S1, S1 e1, class M2, class S2, S2 e2>
constexpr pair<M1, M2> mint_pair() {
return {M1(e1), M2(e2)};
}
#if defined(ATCODER_SEGTREE_HPP)
template <class S, S e>
using segtree_min = segtree<S, op_min<S>, e_const<S, e>>;
template <class S, S e>
using segtree_max = segtree<S, op_max<S>, e_const<S, e>>;
template <class S> using segtree_add = segtree<S, op_add<S>, e_const<S, 0>>;
template <class M, class S = int>
using segtree_add_mint = segtree<M, op_add<M>, mint_const<M, S, 0>>;
#endif
#if defined(ATCODER_LAZYSEGTREE_HPP)
template <class S, S e>
using lazy_segtree_min =
lazy_segtree<S, op_min<S>, e_const<S, e>, pair<S, S>,
mapping_affine<S, S, S>, composition_affine<S, S>,
e_pair<S, 1, S, 0>>;
template <class S, S e>
using lazy_segtree_max =
lazy_segtree<S, op_max<S>, e_const<S, e>, pair<S, S>,
mapping_affine<S, S, S>, composition_affine<S, S>,
e_pair<S, 1, S, 0>>;
template <class S>
using lazy_segtree_add =
lazy_segtree<pair<int, S>, op_add_pair<int, S>, e_pair<int, 0, S, 0>,
pair<S, S>, mapping_len_and_affine<S, S, int, S>,
composition_affine<S, S>, e_pair<S, 1, S, 0>>;
template <class M, class S = int>
using lazy_segtree_add_mint =
lazy_segtree<pair<int, M>, op_add_pair<int, M>,
mint_pair<int, int, 0, M, S, 0>, pair<M, M>,
mapping_len_and_affine<M, M, int, M>, composition_affine<M, M>,
mint_pair<M, S, 1, M, S, 0>>;
#endif
// 一次行列多項式 det(M0 + x M1) を係数列として求める。
// 体上でのみ動作する（除算が必要）。M0, M1 は N×N。
// 多項式は a[0] + a[1]x + ... + a[N]x^N の昇順で返す。
// M1 を掃き出して I にし、det(xI + A) を特性多項式に帰着する。
// M1 が特異でも列に x を掛ける操作を挟むことで次数 1 を保つ。
// 計算量 O(N^3)。
#include <cassert>
#include <utility>
#include <vector>
template <class T>
void hessenberg_reduction(std::vector<std::vector<T>> &matrix) {
const int n = static_cast<int>(matrix.size());
assert(n == 0 || static_cast<int>(matrix[0].size()) == n);
for (int i = 1; i < n; ++i)
assert(static_cast<int>(matrix[i].size()) == n);
for (int r = 0; r < n - 2; ++r) {
int piv = -1;
for (int h = r + 1; h < n; ++h) {
if (matrix[h][r] != T()) {
piv = h;
break;
}
}
if (piv < 0)
continue;
if (piv != r + 1) {
matrix[r + 1].swap(matrix[piv]);
for (int i = 0; i < n; ++i)
std::swap(matrix[i][r + 1], matrix[i][piv]);
}
const T rinv = T(1) / matrix[r + 1][r];
for (int i = r + 2; i < n; ++i) {
const T coef = matrix[i][r] * rinv;
if (coef == T())
continue;
for (int j = 0; j < n; ++j)
matrix[i][j] -= matrix[r + 1][j] * coef;
for (int j = 0; j < n; ++j)
matrix[j][r + 1] += matrix[j][i] * coef;
}
}
}
template <class T>
std::vector<T> characteristic_polynomial(std::vector<std::vector<T>> matrix) {
const int n = static_cast<int>(matrix.size());
assert(n == 0 || static_cast<int>(matrix[0].size()) == n);
for (int i = 1; i < n; ++i)
assert(static_cast<int>(matrix[i].size()) == n);
hessenberg_reduction(matrix);
// p[i] = det(x I_i - matrix[0..i-1][0..i-1])（係数は昇順）
std::vector<std::vector<T>> p(n + 1);
p[0] = {T(1)};
for (int i = 0; i < n; ++i) {
p[i + 1].assign(i + 2, T());
for (int j = 0; j <= i; ++j)
p[i + 1][j + 1] += p[i][j];
for (int j = 0; j <= i; ++j)
p[i + 1][j] -= p[i][j] * matrix[i][i];
T betas = T(1);
for (int j = i - 1; j >= 0; --j) {
betas *= matrix[j + 1][j];
const T hb = (T() - matrix[j][i]) * betas;
for (int k = 0; k <= j; ++k)
p[i + 1][k] += hb * p[j][k];
}
}
return p[n];
}
template <class T>
std::vector<T>
determinant_of_linear_matrix_polynomial(std::vector<std::vector<T>> M0,
std::vector<std::vector<T>> M1) {
const int n = static_cast<int>(M0.size());
assert(static_cast<int>(M1.size()) == n);
if (n == 0)
return {T(1)};
for (int i = 0; i < n; ++i) {
assert(static_cast<int>(M0[i].size()) == n);
assert(static_cast<int>(M1[i].size()) == n);
}
int multiply_by_x = 0; // 「特定の列に x を掛ける」操作の回数
T det_inv = T(1); // 1 / (det A det B)
for (int p = 0; p < n; ++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, T());
// x^2 の項を発生させないため、M1[0..p-1][p] を先に消す（列基本変形）。
for (int row = 0; row < p; ++row) {
const T v = M1[row][p];
if (v == T())
continue;
M1[row][p] = T();
for (int i = 0; i < n; ++i)
M0[i][p] -= v * M0[i][row];
}
// (M0 + x M1) の p 列に x を掛ける（M1 の p 列が 0 のとき swap で実現できる）。
for (int i = 0; i < n; ++i)
std::swap(M0[i][p], M1[i][p]);
--p; // 同じ列をやり直す（高々 n 回）
continue;
}
if (pivot != p) {
M0[pivot].swap(M0[p]);
M1[pivot].swap(M1[p]);
det_inv = T() - det_inv; // *= -1
}
const T v = M1[p][p];
assert(v != T());
det_inv *= v;
const T vinv = T(1) / v;
for (int col = 0; col < n; ++col) {
M0[p][col] *= vinv;
M1[p][col] *= vinv;
}
for (int row = 0; row < n; ++row) {
if (row == p)
continue;
const T coef = M1[row][p];
if (coef == T())
continue;
for (int col = 0; col < n; ++col) {
M0[row][col] -= M0[p][col] * coef;
M1[row][col] -= M1[p][col] * coef;
}
}
}
// M1 = I なので det(x I + M0) を求める（特性多項式 det(x I - (-M0))）。
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j)
M0[i][j] = T() - M0[i][j];
}
std::vector<T> poly = characteristic_polynomial(M0);
for (T &c : poly)
c *= det_inv;
if (multiply_by_x > 0)
poly.erase(poly.begin(), poly.begin() + multiply_by_x);
poly.resize(n + 1, T());
return poly;
}
void solve() {
in_d(int, n);
auto m0 = make_vector({n, n}, mint());
in(m0);
auto m1 = make_vector({n, n}, mint());
in(m1);
out_sep("\n");
out(determinant_of_linear_matrix_polynomial(m0, m1));
}
int main(void) {
cin.tie(nullptr);
ios::sync_with_stdio(false);
cout << fixed << setprecision(15);
cerr << fixed << setprecision(15);
int t = 1;
while (t--) {
solve();
}
return 0;
}