#include #ifdef LOCAL #include #else #define debug(...) void(0) #endif template std::istream& operator>>(std::istream& is, std::vector& v) { for (auto& e : v) { is >> e; } return is; } template std::ostream& operator<<(std::ostream& os, const std::vector& v) { for (std::string_view sep = ""; const auto& e : v) { os << std::exchange(sep, " ") << e; } return os; } template bool chmin(T& x, U&& y) { return y < x and (x = std::forward(y), true); } template bool chmax(T& x, U&& y) { return x < y and (x = std::forward(y), true); } template void mkuni(std::vector& v) { std::ranges::sort(v); auto result = std::ranges::unique(v); v.erase(result.begin(), result.end()); } template int lwb(const std::vector& v, const T& x) { return std::distance(v.begin(), std::ranges::lower_bound(v, x)); } #include template std::vector characteristic_polynomial(std::vector> M) { assert(M.empty() or M.size() == M[0].size()); int n = M.size(); // reduce M to upper Hessenberg form for (int j = 0; j < n - 2; j++) { for (int i = j + 2; i < n; i++) { if (M[i][j] != 0) { std::swap(M[j + 1], M[i]); for (int k = 0; k < n; k++) std::swap(M[k][j + 1], M[k][i]); break; } } if (M[j + 1][j] == 0) continue; auto inv = T(1) / M[j + 1][j]; for (int i = j + 2; i < n; i++) { auto coef = M[i][j] * inv; for (int k = j; k < n; k++) M[i][k] -= coef * M[j + 1][k]; for (int k = 0; k < n; k++) M[k][j + 1] += coef * M[k][i]; } } // compute the characteristic polynomial of upper Hessenberg matrix M std::vector> p(n + 1); p[0] = {T(1)}; for (int i = 0; i < n; i++) { p[i + 1].resize(i + 2); for (int j = 0; j <= i; j++) { p[i + 1][j + 1] += p[i][j]; p[i + 1][j] -= p[i][j] * M[i][i]; } T betas = 1; for (int j = i - 1; j >= 0; j--) { betas *= M[j + 1][j]; T coef = -betas * M[j][i]; for (int k = 0; k <= j; k++) p[i + 1][k] += coef * p[j][k]; } } return p[n]; } template std::vector determinant_polynomial(std::vector> M0, std::vector> M1) { assert(M0.size() == M1.size()); assert(M0.size() == M0[0].size()); assert(M1.size() == M1[0].size()); int n = M0.size(), off = 0; T prod = 1; for (int p = 0; p < n; p++) { int pivot = -1; for (int i = p; i < n; i++) { if (M1[i][p] != 0) { pivot = i; break; } } if (pivot == -1) { if (++off > n) return std::vector(n + 1, 0); for (int i = 0; i < p; i++) { for (int k = 0; k < n; k++) M0[k][p] -= M1[i][p] * M0[k][i]; M1[i][p] = 0; } for (int i = 0; i < n; i++) std::swap(M0[i][p], M1[i][p]); p--; continue; } if (pivot != p) { std::swap(M0[p], M0[pivot]); std::swap(M1[p], M1[pivot]); prod *= -1; } prod *= M1[p][p]; auto inv = T(1) / M1[p][p]; for (int j = 0; j < n; j++) { M0[p][j] *= inv; M1[p][j] *= inv; } for (int i = 0; i < n; i++) { if (i == p) continue; auto coef = M1[i][p]; for (int j = 0; j < n; j++) { M0[i][j] -= M0[p][j] * coef; M1[i][j] -= M1[p][j] * coef; } } } for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { M0[i][j] *= -1; } } auto poly = characteristic_polynomial(M0); std::vector res(n + 1, 0); for (int i = 0; i + off <= n; i++) res[i] = prod * poly[i + off]; return res; } using ll = long long; using namespace std; using mint = atcoder::modint998244353; int main() { ios::sync_with_stdio(false); cin.tie(nullptr); cout << fixed << setprecision(15); int H, W; cin >> H >> W; vector A(H, vector(W)), B(H, vector(W)); cin >> A >> B; vector C(H, vector(H, 0)), D(H, vector(H, 0)); for (int i = 0; i < H; i++) { D[i][i] = 1; for (int j = 0; j < H; j++) { for (int k = 0; k < W; k++) { C[i][j] += mint(A[i][k]) * B[j][k]; } } } auto res = determinant_polynomial(C, D); mint ans = accumulate(res.begin(), res.end(), mint(0)) - 1; cout << ans.val() << "\n"; return 0; }