This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub suisen-cp/cp-library-cpp
#define PROBLEM "https://judge.yosupo.jp/problem/characteristic_polynomial" #include <iostream> #include <atcoder/modint> using mint = atcoder::modint998244353; std::istream& operator>>(std::istream& in, mint &a) { long long e; in >> e; a = e; return in; } std::ostream& operator<<(std::ostream& out, const mint &a) { out << a.val(); return out; } #include "library/linear_algebra/characteristic_polynomial.hpp" using suisen::Matrix; using suisen::characteristic_polynomial; int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); int n; std::cin >> n; Matrix<mint> A(n); for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) std::cin >> A[i][j]; std::vector<mint> p = characteristic_polynomial(A); for (int i = 0; i <= n; ++i) std::cout << p[i] << " \n"[i == n]; return 0; }
#line 1 "test/src/linear_algebra/characteristic_polynomial/characteristic_polynomial.test.cpp" #define PROBLEM "https://judge.yosupo.jp/problem/characteristic_polynomial" #include <iostream> #include <atcoder/modint> using mint = atcoder::modint998244353; std::istream& operator>>(std::istream& in, mint &a) { long long e; in >> e; a = e; return in; } std::ostream& operator<<(std::ostream& out, const mint &a) { out << a.val(); return out; } #line 1 "library/linear_algebra/characteristic_polynomial.hpp" #line 1 "library/linear_algebra/hessenberg_reduction.hpp" #line 1 "library/linear_algebra/matrix.hpp" #include <algorithm> #include <cassert> #include <optional> #include <vector> namespace suisen { template <typename T> struct Matrix { std::vector<std::vector<T>> dat; Matrix() = default; Matrix(int n) : Matrix(n, n) {} Matrix(int n, int m, T fill_value = T(0)) : dat(n, std::vector<T>(m, fill_value)) {} Matrix(const std::vector<std::vector<T>>& dat) : dat(dat) {} const std::vector<T>& operator[](int i) const { return dat[i]; } std::vector<T>& operator[](int i) { return dat[i]; } operator std::vector<std::vector<T>>() const { return dat; } friend bool operator==(const Matrix<T>& A, const Matrix<T>& B) { return A.dat == B.dat; } friend bool operator!=(const Matrix<T>& A, const Matrix<T>& B) { return A.dat != B.dat; } std::pair<int, int> shape() const { return dat.empty() ? std::make_pair<int, int>(0, 0) : std::make_pair<int, int>(dat.size(), dat[0].size()); } int row_size() const { return dat.size(); } int col_size() const { return dat.empty() ? 0 : dat[0].size(); } friend Matrix<T>& operator+=(Matrix<T>& A, const Matrix<T>& B) { assert(A.shape() == B.shape()); auto [n, m] = A.shape(); for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) A.dat[i][j] += B.dat[i][j]; return A; } friend Matrix<T>& operator-=(Matrix<T>& A, const Matrix<T>& B) { assert(A.shape() == B.shape()); auto [n, m] = A.shape(); for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) A.dat[i][j] -= B.dat[i][j]; return A; } friend Matrix<T>& operator*=(Matrix<T>& A, const Matrix<T>& B) { return A = A * B; } friend Matrix<T>& operator*=(Matrix<T>& A, const T& val) { for (auto& row : A.dat) for (auto& elm : row) elm *= val; return A; } friend Matrix<T>& operator/=(Matrix<T>& A, const T& val) { return A *= T(1) / val; } friend Matrix<T>& operator/=(Matrix<T>& A, const Matrix<T>& B) { return A *= B.inv(); } friend Matrix<T> operator+(Matrix<T> A, const Matrix<T>& B) { A += B; return A; } friend Matrix<T> operator-(Matrix<T> A, const Matrix<T>& B) { A -= B; return A; } friend Matrix<T> operator*(const Matrix<T>& A, const Matrix<T>& B) { assert(A.col_size() == B.row_size()); const int n = A.row_size(), m = A.col_size(), l = B.col_size(); if (std::min({ n, m, l }) <= 70) { // naive Matrix<T> C(n, l); for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) for (int k = 0; k < l; ++k) { C.dat[i][k] += A.dat[i][j] * B.dat[j][k]; } return C; } // strassen const int nl = 0, nm = n >> 1, nr = nm + nm; const int ml = 0, mm = m >> 1, mr = mm + mm; const int ll = 0, lm = l >> 1, lr = lm + lm; auto A00 = A.submatrix(nl, nm, ml, mm), A01 = A.submatrix(nl, nm, mm, mr); auto A10 = A.submatrix(nm, nr, ml, mm), A11 = A.submatrix(nm, nr, mm, mr); auto B00 = B.submatrix(ml, mm, ll, lm), B01 = B.submatrix(ml, mm, lm, lr); auto B10 = B.submatrix(mm, mr, ll, lm), B11 = B.submatrix(mm, mr, lm, lr); auto P0 = (A00 + A11) * (B00 + B11); auto P1 = (A10 + A11) * B00; auto P2 = A00 * (B01 - B11); auto P3 = A11 * (B10 - B00); auto P4 = (A00 + A01) * B11; auto P5 = (A10 - A00) * (B00 + B01); auto P6 = (A01 - A11) * (B10 + B11); Matrix<T> C(n, l); C.assign_submatrix(nl, ll, P0 + P3 - P4 + P6), C.assign_submatrix(nl, lm, P2 + P4); C.assign_submatrix(nm, ll, P1 + P3), C.assign_submatrix(nm, lm, P0 + P2 - P1 + P5); // fractions if (l != lr) { for (int i = 0; i < nr; ++i) for (int j = 0; j < mr; ++j) C.dat[i][lr] += A.dat[i][j] * B.dat[j][lr]; } if (m != mr) { for (int i = 0; i < nr; ++i) for (int k = 0; k < l; ++k) C.dat[i][k] += A.dat[i][mr] * B.dat[mr][k]; } if (n != nr) { for (int j = 0; j < m; ++j) for (int k = 0; k < l; ++k) C.dat[nr][k] += A.dat[nr][j] * B.dat[j][k]; } return C; } friend Matrix<T> operator*(const T& val, Matrix<T> A) { A *= val; return A; } friend Matrix<T> operator*(Matrix<T> A, const T& val) { A *= val; return A; } friend Matrix<T> operator/(const Matrix<T>& A, const Matrix<T>& B) { return A * B.inv(); } friend Matrix<T> operator/(Matrix<T> A, const T& val) { A /= val; return A; } friend Matrix<T> operator/(const T& val, const Matrix<T>& A) { return val * A.inv(); } friend std::vector<T> operator*(const Matrix<T>& A, const std::vector<T>& x) { const auto [n, m] = A.shape(); assert(m == int(x.size())); std::vector<T> b(n, T(0)); for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) b[i] += A.dat[i][j] * x[j]; return b; } static Matrix<T> e0(int n) { return Matrix<T>(n, Identity::ADD); } static Matrix<T> e1(int n) { return Matrix<T>(n, Identity::MUL); } Matrix<T> pow(long long b) const { assert_square(); assert(b >= 0); Matrix<T> res = e1(row_size()), p = *this; for (; b; b >>= 1) { if (b & 1) res *= p; p *= p; } return res; } Matrix<T> inv() const { return *safe_inv(); } std::optional<Matrix<T>> safe_inv() const { assert_square(); Matrix<T> A = *this; const int n = A.row_size(); for (int i = 0; i < n; ++i) { A[i].resize(2 * n, T{ 0 }); A[i][n + i] = T{ 1 }; } for (int i = 0; i < n; ++i) { for (int k = i; k < n; ++k) if (A[k][i] != T{ 0 }) { std::swap(A[i], A[k]); T c = T{ 1 } / A[i][i]; for (int j = i; j < 2 * n; ++j) A[i][j] *= c; break; } if (A[i][i] == T{ 0 }) return std::nullopt; for (int k = 0; k < n; ++k) if (k != i and A[k][i] != T{ 0 }) { T c = A[k][i]; for (int j = i; j < 2 * n; ++j) A[k][j] -= c * A[i][j]; } } for (auto& row : A.dat) row.erase(row.begin(), row.begin() + n); return A; } T det() const { assert_square(); Matrix<T> A = *this; bool sgn = false; const int n = A.row_size(); for (int j = 0; j < n; ++j) for (int i = j + 1; i < n; ++i) if (A[i][j] != T{ 0 }) { std::swap(A[j], A[i]); T q = A[i][j] / A[j][j]; for (int k = j; k < n; ++k) A[i][k] -= A[j][k] * q; sgn = not sgn; } T res = sgn ? T{ -1 } : T{ +1 }; for (int i = 0; i < n; ++i) res *= A[i][i]; return res; } T det_arbitrary_mod() const { assert_square(); Matrix<T> A = *this; bool sgn = false; const int n = A.row_size(); for (int j = 0; j < n; ++j) for (int i = j + 1; i < n; ++i) { for (; A[i][j].val(); sgn = not sgn) { std::swap(A[j], A[i]); T q = A[i][j].val() / A[j][j].val(); for (int k = j; k < n; ++k) A[i][k] -= A[j][k] * q; } } T res = sgn ? -1 : +1; for (int i = 0; i < n; ++i) res *= A[i][i]; return res; } void assert_square() const { assert(row_size() == col_size()); } Matrix<T> submatrix(int row_begin, int row_end, int col_begin, int col_end) const { Matrix<T> A(row_end - row_begin, col_end - col_begin); for (int i = row_begin; i < row_end; ++i) for (int j = col_begin; j < col_end; ++j) { A[i - row_begin][j - col_begin] = dat[i][j]; } return A; } void assign_submatrix(int row_begin, int col_begin, const Matrix<T>& A) { const int n = A.row_size(), m = A.col_size(); assert(row_begin + n <= row_size() and col_begin + m <= col_size()); for (int i = 0; i < n; ++i) for (int j = 0; j < m; ++j) { dat[row_begin + i][col_begin + j] = A[i][j]; } } private: enum class Identity { ADD, MUL }; Matrix(int n, Identity ident) : Matrix<T>::Matrix(n) { if (ident == Identity::MUL) for (int i = 0; i < n; ++i) dat[i][i] = 1; } }; } // namespace suisen #line 5 "library/linear_algebra/hessenberg_reduction.hpp" namespace suisen { /** * Reference: http://www.phys.uri.edu/nigh/NumRec/bookfpdf/f11-5.pdf * returns H := P^(-1) A P, where H is hessenberg matrix */ template <typename T> Matrix<T> hessenberg_reduction(Matrix<T> A) { A.assert_square(); const int n = A.row_size(); for (int r = 0; r < n - 2; ++r) { int pivot = -1; for (int r2 = r + 1; r2 < n; ++r2) if (A[r2][r] != 0) { pivot = r2; break; } if (pivot < 0) continue; if (pivot != r + 1) { for (int k = 0; k < n; ++k) std::swap(A[r + 1][k], A[pivot][k]); for (int k = 0; k < n; ++k) std::swap(A[k][r + 1], A[k][pivot]); } const T den = T{1} / A[r + 1][r]; for (int r2 = r + 2; r2 < n; ++r2) if (T coef = A[r2][r] * den; coef != 0) { for (int k = r; k < n; ++k) A[r2][k] -= coef * A[r + 1][k]; for (int k = 0; k < n; ++k) A[k][r + 1] += coef * A[k][r2]; } } return A; } } // namespace suisen #line 5 "library/linear_algebra/characteristic_polynomial.hpp" namespace suisen { /** * Reference: https://ipsen.math.ncsu.edu/ps/charpoly3.pdf * returns p(λ) = det(λE - A) */ template <typename T> std::vector<T> characteristic_polynomial(const Matrix<T> &A) { A.assert_square(); const int n = A.row_size(); if (n == 0) return { T{1} }; auto H = hessenberg_reduction(A); /** * +- -+ * | a0 * * * * | * | b1 a1 * * * | * H = | 0 b2 a2 * * | * | 0 0 b3 a3 * | * | 0 0 0 b4 a4 | * +- -+ * p_i(λ) := det(λ*E_i - H[:i][:i]) * p_0(λ) = 1, * p_1(λ) = λ-a_0, * p_i(λ) = (λ-a_{i-1}) p_{i-1}(λ) - Σ[j=0,i-1] p_j(λ) * H_{j,i} * Π[k=j+1,i] b_k. */ std::vector<std::vector<T>> p(n + 1); p[0] = { T{1} }, p[1] = { { -H[0][0], T{1} } }; for (int i = 1; i < n; ++i) { p[i + 1].resize(i + 2, T{0}); for (int k = 0; k < i + 1; ++k) { p[i + 1][k] -= H[i][i] * p[i][k]; p[i + 1][k + 1] += p[i][k]; } T prod_b = T{1}; for (int j = i - 1; j >= 0; --j) { prod_b *= H[j + 1][j]; T coef = H[j][i] * prod_b; for (int k = 0; k < j + 1; ++k) p[i + 1][k] -= coef * p[j][k]; } } return p[n]; } } // namespace suisen #line 19 "test/src/linear_algebra/characteristic_polynomial/characteristic_polynomial.test.cpp" using suisen::Matrix; using suisen::characteristic_polynomial; int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); int n; std::cin >> n; Matrix<mint> A(n); for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) std::cin >> A[i][j]; std::vector<mint> p = characteristic_polynomial(A); for (int i = 0; i <= n; ++i) std::cout << p[i] << " \n"[i == n]; return 0; }