#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define int ll #define INT128_MAX (__int128)(((unsigned __int128) 1 << ((sizeof(__int128) * __CHAR_BIT__) - 1)) - 1) #define INT128_MIN (-INT128_MAX - 1) #define pb push_back #define eb emplace_back #define clock chrono::steady_clock::now().time_since_epoch().count() using namespace std; template ostream& print_tuple(ostream& os, const tuple tu) { os << get(tu); if constexpr (I + 1 != sizeof...(args)) { os << ' '; print_tuple(os, tu); } return os; } template ostream& operator<<(ostream& os, const tuple tu) { return print_tuple(os, tu); } template ostream& operator<<(ostream& os, const pair pr) { return os << pr.first << ' ' << pr.second; } template ostream& operator<<(ostream& os, const array &arr) { for(size_t i = 0; T x : arr) { os << x; if (++i != N) os << ' '; } return os; } template ostream& operator<<(ostream& os, const vector &vec) { for(size_t i = 0; T x : vec) { os << x; if (++i != size(vec)) os << ' '; } return os; } template ostream& operator<<(ostream& os, const set &s) { for(size_t i = 0; T x : s) { os << x; if (++i != size(s)) os << ' '; } return os; } template ostream& operator<<(ostream& os, const multiset &s) { for(size_t i = 0; T x : s) { os << x; if (++i != size(s)) os << ' '; } return os; } template ostream& operator<<(ostream& os, const map &m) { for(size_t i = 0; pair x : m) { os << x.first << " : " << x.second; if (++i != size(m)) os << ", "; } return os; } #ifdef DEBUG #define dbg(...) cerr << '(', _do(#__VA_ARGS__), cerr << ") = ", _do2(__VA_ARGS__) template void _do(T &&x) { cerr << x; } template void _do(T &&x, S&&...y) { cerr << x << ", "; _do(y...); } template void _do2(T &&x) { cerr << x << endl; } template void _do2(T &&x, S&&...y) { cerr << x << ", "; _do2(y...); } #else #define dbg(...) #endif using ll = long long; using ull = unsigned long long; using ldb = long double; using pii = pair; using pll = pair; //#define double ldb template using vc = vector; template using vvc = vc>; template using vvvc = vc>; using vi = vc; using vll = vc; using vvi = vvc; using vvll = vvc; template using min_heap = priority_queue, greater>; template using max_heap = priority_queue; template concept R_invocable = requires(F&& f, Args&&... args) { { std::invoke(std::forward(f), std::forward(args)...) } -> std::same_as; }; template, typename F> requires R_invocable void pSum(rng &&v, F f) { if (!v.empty()) for(T p = *v.begin(); T &x : v | views::drop(1)) x = p = f(p, x); } template> void pSum(rng &&v) { if (!v.empty()) for(T p = *v.begin(); T &x : v | views::drop(1)) x = p = p + x; } template void Unique(rng &v) { ranges::sort(v); v.resize(unique(v.begin(), v.end()) - v.begin()); } template rng invPerm(rng p) { rng ret = p; for(int i = 0; i < ssize(p); i++) ret[p[i]] = i; return ret; } template vi argSort(rng p) { vi id(size(p)); iota(id.begin(), id.end(), 0); ranges::sort(id, {}, [&](int i) { return pair(p[i], i); }); return id; } template, typename F> requires invocable vi argSort(rng p, F proj) { vi id(size(p)); iota(id.begin(), id.end(), 0); ranges::sort(id, {}, [&](int i) { return pair(proj(p[i]), i); }); return id; } template vvi read_graph(int n, int m, int base) { vvi g(n); for(int i = 0; i < m; i++) { int u, v; cin >> u >> v; u -= base, v -= base; g[u].emplace_back(v); if constexpr (!directed) g[v].emplace_back(u); } return g; } template vvi adjacency_list(int n, vc e, int base) { vvi g(n); for(auto [u, v] : e) { u -= base, v -= base; g[u].emplace_back(v); if constexpr (!directed) g[v].emplace_back(u); } return g; } template vc equal_subarrays(vc &v) { vc lr; for(int i = 0, j = 0; i < ssize(v); i = j) { while(j < ssize(v) and v[i] == v[j]) j++; lr.eb(i, j); } return lr; } template requires invocable vc equal_subarrays(vc &v, F proj) { vc lr; for(int i = 0, j = 0; i < ssize(v); i = j) { while(j < ssize(v) and proj(v[i]) == proj(v[j])) j++; lr.eb(i, j); } return lr; } template void setBit(T &msk, int bit, bool x) { (msk &= ~(T(1) << bit)) |= T(x) << bit; } template void onBit(T &msk, int bit) { setBit(msk, bit, true); } template void offBit(T &msk, int bit) { setBit(msk, bit, false); } template void flipBit(T &msk, int bit) { msk ^= T(1) << bit; } template bool getBit(T msk, int bit) { return msk >> bit & T(1); } template T floorDiv(T a, T b) { if (b < 0) a *= -1, b *= -1; return a >= 0 ? a / b : (a - b + 1) / b; } template T ceilDiv(T a, T b) { if (b < 0) a *= -1, b *= -1; return a >= 0 ? (a + b - 1) / b : a / b; } template bool chmin(T &a, T b) { return a > b ? a = b, 1 : 0; } template bool chmax(T &a, T b) { return a < b ? a = b, 1 : 0; } //reference: https://github.com/NyaanNyaan/library/blob/master/modint/montgomery-modint.hpp#L10 //note: mod should be an odd prime less than 2^30. template struct MontgomeryModInt { using mint = MontgomeryModInt; using i32 = int32_t; using u32 = uint32_t; using u64 = uint64_t; static constexpr u32 get_r() { u32 res = 1, base = mod; for(i32 i = 0; i < 31; i++) res *= base, base *= base; return -res; } static constexpr u32 get_mod() { return mod; } static constexpr u32 n2 = -u64(mod) % mod; //2^64 % mod static constexpr u32 r = get_r(); //-P^{-1} % 2^32 u32 a; static u32 reduce(const u64 &b) { return (b + u64(u32(b) * r) * mod) >> 32; } static u32 transform(const u64 &b) { return reduce(u64(b) * n2); } MontgomeryModInt() : a(0) {} MontgomeryModInt(const int64_t &b) : a(transform(b % mod + mod)) {} mint pow(u64 k) const { mint res(1), base(*this); while(k) { if (k & 1) res *= base; base *= base, k >>= 1; } return res; } mint inverse() const { return (*this).pow(mod - 2); } u32 get() const { u32 res = reduce(a); return res >= mod ? res - mod : res; } mint& operator+=(const mint &b) { if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod; return *this; } mint& operator-=(const mint &b) { if (i32(a -= b.a) < 0) a += 2 * mod; return *this; } mint& operator*=(const mint &b) { a = reduce(u64(a) * b.a); return *this; } mint& operator/=(const mint &b) { a = reduce(u64(a) * b.inverse().a); return *this; } mint operator-() { return mint() - mint(*this); } bool operator==(mint b) const { return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a); } bool operator!=(mint b) const { return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a); } friend mint operator+(mint c, mint d) { return c += d; } friend mint operator-(mint c, mint d) { return c -= d; } friend mint operator*(mint c, mint d) { return c *= d; } friend mint operator/(mint c, mint d) { return c /= d; } friend ostream& operator<<(ostream& os, const mint& b) { return os << b.get(); } friend istream& operator>>(istream& is, mint& b) { int64_t val; is >> val; b = mint(val); return is; } }; //using mint = MontgomeryModInt<1'000'000'007>; using mint = MontgomeryModInt<998'244'353>; //#include template struct binomial { vector _fac, _facInv; binomial(int size) : _fac(size), _facInv(size) { assert(size <= (int)Mint::get_mod()); _fac[0] = 1; for(int i = 1; i < size; i++) _fac[i] = _fac[i - 1] * i; if (size > 0) _facInv.back() = 1 / _fac.back(); for(int i = size - 2; i >= 0; i--) _facInv[i] = _facInv[i + 1] * (i + 1); } Mint fac(int i) { return i < 0 ? 0 : _fac[i]; } Mint faci(int i) { return i < 0 ? 0 : _facInv[i]; } Mint inv(int i) { return _facInv[i] * _fac[i - 1]; } Mint binom(int n, int r) { return r < 0 or n < r ? 0 : fac(n) * faci(r) * faci(n - r); } Mint catalan(int i) { return binom(2 * i, i) - binom(2 * i, i + 1); } Mint excatalan(int n, int m, int k) { //(+1) * n, (-1) * m, prefix sum > -k if (k > m) return binom(n + m, m); else if (k > m - n) return binom(n + m, m) - binom(n + m, m - k); else return Mint(0); } }; template struct matrix : vector> { matrix(int n, int m) : vector>(n, vector(m, 0)) {} matrix(int n) : vector>(n, vector(n, 0)) {} matrix(vvc M) : vvc(M) {} int n() const { return ssize(*this); } int m() const { return n() == 0 ? 0 : ssize((*this)[0]); } static matrix I(int n) { auto res = matrix(n, n); for(int i = 0; i < n; i++) res[i][i] = 1; return res; } matrix& operator+=(const matrix &b) { assert(n() == b.n()); assert(m() == b.m()); for(int i = 0; i < n(); i++) for(int j = 0; j < m(); j++) (*this)[i][j] += b[i][j]; return *this; } matrix& operator-=(const matrix &b) { assert(n() == b.n()); assert(m() == b.m()); for(int i = 0; i < n(); i++) for(int j = 0; j < m(); j++) (*this)[i][j] -= b[i][j]; return *this; } matrix& operator*=(const matrix &b) { assert(m() == b.n()); auto res = matrix(n(), b.m()); for(int i = 0; i < n(); i++) for(int k = 0; k < m(); k++) for(int j = 0; j < b.m(); j++) res[i][j] += (*this)[i][k] * b[k][j]; this -> swap(res); return *this; } matrix pow(ll k) const { assert(n() == m()); auto res = I(n()), base = *this; while(k) { if (k & 1) res *= base; base *= base, k >>= 1; } return res; } tuple, int> eliminate() { int sgn = 1; matrix M = *this; vector pivot_row; for(int row = 0, col = 0; row < n() and col < m(); col++) { int p_row = -1; for(int i = row; i < n() and p_row == -1; i++) if (M[i][col] != 0) p_row = i; if (p_row == -1) continue; pivot_row.emplace_back(row); if (row != p_row) { for(int j = col; j < m(); j++) swap(M[row][j], M[p_row][j]); sgn *= -1; } for(int i = 0; i < n(); i++) { if (i == row or M[i][col] == 0) continue; Mint s = M[i][col] / M[row][col]; for(int j = col; j < m(); j++) M[i][j] -= M[row][j] * s; } row++; } return {M, pivot_row, sgn}; } Mint det() { assert(n() == m()); auto [M, pr, sgn] = eliminate(); if (ssize(pr) != n()) { return Mint(0); } else { Mint d = sgn; for(int i = 0; i < n(); i++) d *= M[i][i]; return d; } } int rank() { return get<1>(eliminate()).size(); } pair inv() { assert(n() == m()); matrix M(n(), 2 * n()); for(int i = 0; i < n(); i++) { for(int j = 0; j < n(); j++) M[i][j] = (*this)[i][j]; M[i][n() + i] = 1; } matrix tmp = get<0>(M.eliminate()); matrix MI(n(), n()); for(int i = 0; i < n(); i++) { if (tmp[i][i] == 0) return {false, matrix(0, 0)}; Mint r = tmp[i][i].inverse(); for(int j = 0; j < n(); j++) MI[i][j] = tmp[i][j + n()] * r; } return {true, MI}; } pair, matrix> solve_linear(vector b) { assert(n() == ssize(b)); matrix M(n(), m() + 1); for(int i = 0; i < n(); i++) { for(int j = 0; j < m(); j++) M[i][j] = (*this)[i][j]; M[i][m()] = b[i]; } auto [N, pr, _] = M.eliminate(); vector x(m()); vector where(m(), -1), inv_where(m(), -1); for(int row : pr) { int col = 0; while(N[row][col] == 0) col++; if (col < m()) where[col] = row, inv_where[row] = col; } for(int i = 0; i < m(); i++) if (where[i] != -1) x[i] = N[where[i]][m()] / N[where[i]][i]; for(int i = 0; i < n(); i++) { Mint s = -N[i][m()]; for(int j = 0; j < m(); j++) s += x[j] * N[i][j]; if (s != Mint(0)) return {vector(), matrix(0)}; } matrix basis(m() - ssize(pr), m()); for(int col = 0, last_row = 0, k = 0; col < m(); col++) { if (where[col] != -1) { last_row = where[col]; } else { basis[k][col] = 1; for(int i = 0; i <= last_row; i++) basis[k][inv_where[i]] = -N[i][col] / N[i][inv_where[i]]; k++; } } return {x, basis}; } matrix operator-() { return matrix(n(), m()) - (*this); } friend matrix operator+(matrix a, matrix b) { return a += b; } friend matrix operator-(matrix a, matrix b) { return a -= b; } friend matrix operator*(matrix a, matrix b) { return a *= b; } friend ostream& operator<<(ostream& os, const matrix& b) { for(int i = 0; i < b.n(); i++) { os << '\n'; for(int j = 0; j < b.m(); j++) os << b[i][j] << ' '; } return os; } friend istream& operator>>(istream& is, matrix& b) { for(int i = 0; i < b.n(); i++) for(int j = 0; j < b.m(); j++) is >> b[i][j]; return is; } }; template array operator+(array a, array b) { return {a[0] + b[0], a[1] + b[1]}; } template array operator-(array a, array b) { return {a[0] - b[0], a[1] - b[1]}; } template array operator*(Mint m, array b) { return {m * b[0], m * b[1]}; } template vc characteristic_polynomial(vvc M) { if (M.empty()) return {1}; assert(size(M) == size(M[0])); const int N = size(M); vc P(N, vc>(N)); for(int i = 0; i < N; i++) for(int j = 0; j < N; j++) P[i][j] = {-M[i][j], i == j}; //reduce to lower Hessenberg Matrix for(int r = 0; r < N - 2; r++) { const int c = r + 1; int i = c; while(i < N and P[r][i][0] == 0) i++; if (i == N) continue; if (i > c) { P[i].swap(P[c]); for(int j = 0; j < N; j++) swap(P[i][j], P[c][j]); } Mint inv = P[r][c][0].inverse(); for(int j = c + 1; j < N; j++) { Mint R = -inv * P[r][j][0]; for(int k = 0; k < N; k++) P[k][j] = P[k][j] + R * P[k][c]; } for(int j = r + 2; j < N; j++) { Mint R = -P[r + 1][j][1]; for(int k = 0; k < N; k++) P[r + 1][k] = P[r + 1][k] + R * P[j][k]; } } auto add = [&](vector a, vector b) { if (ssize(a) < ssize(b)) a.resize(size(b)); for(int i = 0; i < ssize(b); i++) a[i] += b[i]; return a; }; auto mul = [&](vector &a, array b) { vector c(ssize(a) + 1); for(int i = 0; i < ssize(a); i++) { c[i] += a[i] * b[0]; c[i + 1] += a[i] * b[1]; } return c; }; auto neg = [&](array a) { return array{-a[0], -a[1]}; }; //DP vc dp(1, vc{1}); for(int i = 0; i < N - 1; i++) { vc> nxt(i + 2); for(int j = 0; j < ssize(dp); j++) { nxt[j] = add(nxt[j], mul(dp[j], P[i][i + 1])); nxt[i + 1] = add(nxt[i + 1], mul(dp[j], (i - j) % 2 ? neg(P[i][j]) : P[i][j])); } dp.swap(nxt); } vc C = {0}; for(int j = 0; j < N; j++) C = add(C, mul(dp[j], (N - 1 - j) % 2 ? neg(P[N - 1][j]) : P[N - 1][j])); return C; } //compute det(M0 + M1x) with N/MOD probability to fail template auto det_poly(matrix M0, matrix M1) { const int N = ssize(M0); assert(N > 0 and ssize(M1) == N and ssize(M0[0]) == N and ssize(M1[0]) == N); mt19937_64 rng(clock); mint a = rng() % Mint::get_mod(); for(int i = 0; i < N; i++) for(int j = 0; j < N; j++) M0[i][j] += a * M1[i][j]; swap(M0, M1); auto [ok, M1I] = M1.inv(); if (!ok) return pair(false, vc(N + 1, 0)); auto p = characteristic_polynomial(-M1I * M0); mint inv_det = M1I.det().inverse(); for(mint &x : p) x *= inv_det; vc q(N + 1); for(int i = 0; i < ssize(p); i++) q[N - i] = p[i]; binomial bn(N + 1); vc r(N + 1); for(int i = 0; i < ssize(q); i++) { mint mul = 1; for(int j = 0; j <= i; j++, mul *= -a) r[i - j] += q[i] * mul * bn.binom(i, j); } return pair(true, r); } signed main() { ios::sync_with_stdio(false), cin.tie(NULL); int n; cin >> n; matrix M0(n, n), M1(n, n); cin >> M0 >> M1; for(mint x : det_poly(M0, M1).second) cout << x << '\n'; return 0; }