sotanishy's code snippets for competitive programming
View the Project on GitHub sotanishy/cp-library-cpp
#include "math/linalg/characteristic_polynomial.hpp"
正方行列の特性多項式を求める.
Polynomial characteristic_polynomial(SquareMatrix mat)
#pragma once #include <algorithm> #include <vector> #include "../polynomial.cpp" #include "matrix.hpp" template <typename mint> Polynomial<mint> characteristic_polynomial(Matrix<mint> mat) { mat.assert_square(); const int n = mat.row_size(); if (n == 0) return {1}; // stage 1: reduce mat to upper Hessenberg form for (int j = 0; j < n; ++j) { int i = j + 1; while (i < n && mat[i][j] == 0) ++i; if (i == n) continue; if (i != j + 1) { // swap mat[i], mat[j+1] mat[i].swap(mat[j + 1]); // swap mat[:,i], mat[:,j+1] for (int k = 0; k < n; ++k) { std::swap(mat[k][i], mat[k][j + 1]); } } mint inv = mint(1) / mat[j + 1][j]; for (int k = j + 2; k < n; ++k) { mint v = mat[k][j] * inv; // mat[k] -= mat[j+1] * v for (int l = j; l < n; ++l) { mat[k][l] -= mat[j + 1][l] * v; } // mat[:,j+1] += mat[:,k] * v for (int l = 0; l < n; ++l) { mat[l][j + 1] += mat[l][k] * v; } } } // stage 2: recursively compute char polys of leading principal submatrices std::vector<Polynomial<mint>> p(n + 1); p[0] = {1}; for (int i = 0; i < n; ++i) { p[i + 1].resize(i + 1 + 1); // (x - mat[i,i]) * p[i] for (int j = 0; j <= i; ++j) { p[i + 1][j] -= mat[i][i] * p[i][j]; p[i + 1][j + 1] += p[i][j]; } mint beta = 1; for (int k = i - 1; k >= 0; --k) { beta *= mat[k + 1][k]; p[i + 1] -= p[k] * mat[k][i] * beta; } } return p[n]; }
#line 2 "math/linalg/characteristic_polynomial.hpp" #include <algorithm> #include <vector> #line 3 "math/polynomial.cpp" #include <cassert> #line 5 "math/polynomial.cpp" #line 2 "convolution/ntt.hpp" #include <bit> #line 4 "convolution/ntt.hpp" constexpr int get_primitive_root(int mod) { if (mod == 167772161) return 3; if (mod == 469762049) return 3; if (mod == 754974721) return 11; if (mod == 998244353) return 3; if (mod == 1224736769) return 3; } template <typename mint> void ntt(std::vector<mint>& a) { constexpr int mod = mint::mod(); constexpr mint primitive_root = get_primitive_root(mod); const int n = a.size(); for (int m = n; m > 1; m >>= 1) { mint omega = primitive_root.pow((mod - 1) / m); for (int s = 0; s < n / m; ++s) { mint w = 1; for (int i = 0; i < m / 2; ++i) { mint l = a[s * m + i]; mint r = a[s * m + i + m / 2]; a[s * m + i] = l + r; a[s * m + i + m / 2] = (l - r) * w; w *= omega; } } } } template <typename mint> void intt(std::vector<mint>& a) { constexpr int mod = mint::mod(); constexpr mint primitive_root = get_primitive_root(mod); const int n = a.size(); for (int m = 2; m <= n; m <<= 1) { mint omega = primitive_root.pow((mod - 1) / m).inv(); for (int s = 0; s < n / m; ++s) { mint w = 1; for (int i = 0; i < m / 2; ++i) { mint l = a[s * m + i]; mint r = a[s * m + i + m / 2] * w; a[s * m + i] = l + r; a[s * m + i + m / 2] = l - r; w *= omega; } } } } template <typename mint> std::vector<mint> convolution(std::vector<mint> a, std::vector<mint> b) { const int size = a.size() + b.size() - 1; const int n = std::bit_ceil((unsigned int)size); a.resize(n); b.resize(n); ntt(a); ntt(b); for (int i = 0; i < n; ++i) a[i] *= b[i]; intt(a); a.resize(size); mint n_inv = mint(n).inv(); for (int i = 0; i < size; ++i) a[i] *= n_inv; return a; } #line 7 "math/polynomial.cpp" template <typename mint> class Polynomial : public std::vector<mint> { using Poly = Polynomial; public: using std::vector<mint>::vector; using std::vector<mint>::operator=; Poly pre(int size) const { return Poly(this->begin(), this->begin() + std::min((int)this->size(), size)); } Poly rev(int deg = -1) const { auto ret = *this; if (deg != -1) ret.resize(deg, 0); return Poly(ret.rbegin(), ret.rend()); } void trim() { while (!this->empty() && this->back() == 0) this->pop_back(); } // --- unary operation --- Poly& operator-() const { auto ret = *this; for (auto& x : ret) x = -x; return ret; } // -- binary operation with scalar --- Poly& operator+=(const mint& rhs) { if (this->empty()) this->resize(1); (*this)[0] += rhs; return *this; } Poly& operator-=(const mint& rhs) { if (this->empty()) this->resize(1); (*this)[0] -= rhs; return *this; } Poly& operator*=(const mint& rhs) { for (auto& x : *this) x *= rhs; return *this; } Poly& operator/=(const mint& rhs) { return *this *= rhs.inv(); } Poly operator+(const mint& rhs) const { return Poly(*this) += rhs; } Poly operator-(const mint& rhs) const { return Poly(*this) -= rhs; } Poly operator*(const mint& rhs) const { return Poly(*this) *= rhs; } Poly operator/(const mint& rhs) const { return Poly(*this) /= rhs; } // --- binary operation with polynomial --- Poly& operator+=(const Poly& rhs) { if (this->size() < rhs.size()) this->resize(rhs.size()); for (int i = 0; i < (int)rhs.size(); ++i) (*this)[i] += rhs[i]; return *this; } Poly& operator-=(const Poly& rhs) { if (this->size() < rhs.size()) this->resize(rhs.size()); for (int i = 0; i < (int)rhs.size(); ++i) (*this)[i] -= rhs[i]; return *this; } Poly& operator*=(const Poly& rhs) { *this = convolution(*this, rhs); return *this; // // naive convolution O(N^2) // std::vector<mint> res(this->size() + rhs.size() - 1); // for (int i = 0; i < (int)this->size(); ++i) { // for (int j = 0; j < (int)rhs.size(); ++j) { // res[i + j] += (*this)[i] * rhs[j]; // } // } // return *this = res; } Poly& operator/=(const Poly& rhs) { if (this->size() < rhs.size()) { this->clear(); return *this; } int n = this->size() - rhs.size() + 1; return *this = (rev().pre(n) * rhs.rev().inv(n)).pre(n).rev(n); } Poly& operator%=(const Poly& rhs) { *this -= *this / rhs * rhs; trim(); return *this; } std::pair<Poly, Poly> divmod(const Poly& rhs) { auto q = *this / rhs; auto r = *this - q * rhs; r.trim(); return {q, r}; } Poly operator+(const Poly& rhs) const { return Poly(*this) += rhs; } Poly operator-(const Poly& rhs) const { return Poly(*this) -= rhs; } Poly operator*(const Poly& rhs) const { return Poly(*this) *= rhs; } Poly operator/(const Poly& rhs) const { return Poly(*this) /= rhs; } Poly operator%(const Poly& rhs) const { return Poly(*this) %= rhs; } // --- shift operation --- Poly operator<<(int n) const { auto ret = *this; ret.insert(ret.begin(), n, 0); return ret; } Poly operator>>(int n) const { if ((int)this->size() <= n) return {}; auto ret = *this; ret.erase(ret.begin(), ret.begin() + n); return ret; } // --- evaluation --- mint operator()(const mint& x) { mint y = 0, powx = 1; for (int i = 0; i < (int)this->size(); ++i) { for (auto c : *this) { y += c * powx; powx *= x; } return y; } } // --- other operations --- Poly inv(int deg = -1) const { assert((*this)[0] != mint(0)); if (deg == -1) deg = this->size(); Poly res = {(*this)[0].inv()}; for (int d = 1; d < deg; d <<= 1) { auto f = pre(2 * d); auto g = res; f.resize(2 * d); g.resize(2 * d); // g_{n+1} = g_n * (2 - g_n * f) mod x^{2^{n+1}} ntt(f); ntt(g); for (int i = 0; i < 2 * d; ++i) f[i] *= g[i]; intt(f); for (int i = 0; i < d; ++i) f[i] = 0; ntt(f); for (int i = 0; i < 2 * d; ++i) f[i] *= g[i]; intt(f); res.resize(2 * d); auto coef = mint(2 * d).inv().pow(2); for (int i = d; i < 2 * d; ++i) res[i] = -f[i] * coef; } return res.pre(deg); } Poly exp(int deg = -1) const { assert((*this)[0] == mint(0)); if (deg == -1) deg = this->size(); Poly ret = {mint(1)}; for (int i = 1; i < deg; i <<= 1) { ret = (ret * (this->pre(i << 1) + mint(1) - ret.log(i << 1))) .pre(i << 1); } return ret; } Poly log(int deg = -1) const { assert((*this)[0] == mint(1)); if (deg == -1) deg = this->size(); return (diff() * inv(deg)).pre(deg - 1).integral(); } Poly pow(long long k, int deg = -1) const { if (k == 0) return {1}; if (deg == -1) deg = this->size(); auto ret = *this; int cnt0 = 0; while (cnt0 < (int)ret.size() && ret[cnt0] == 0) ++cnt0; if (cnt0 > (deg - 1) / k) return {}; ret = ret >> cnt0; deg -= cnt0 * k; ret = ((ret / ret[0]).log(deg) * k).exp(deg) * ret[0].pow(k); ret = ret << (cnt0 * k); return ret; } Poly diff() const { Poly ret(std::max(0, (int)this->size() - 1)); for (int i = 1; i <= (int)ret.size(); ++i) ret[i - 1] = (*this)[i] * mint(i); return ret; } Poly integral() const { Poly ret(this->size() + 1); ret[0] = mint(0); for (int i = 0; i < (int)ret.size() - 1; ++i) ret[i + 1] = (*this)[i] / mint(i + 1); return ret; } Poly taylor_shift(long long c) const { const int n = this->size(); std::vector<mint> fact(n, 1), fact_inv(n, 1); for (int i = 1; i < n; ++i) fact[i] = fact[i - 1] * i; fact_inv[n - 1] = mint(1) / fact[n - 1]; for (int i = n - 1; i > 0; --i) fact_inv[i - 1] = fact_inv[i] * i; auto ret = *this; Poly e(n + 1); e[0] = 1; mint p = c; for (int i = 1; i < n; ++i) { ret[i] *= fact[i]; e[i] = p * fact_inv[i]; p *= c; } ret = (ret.rev() * e).pre(n).rev(); for (int i = n - 1; i >= 0; --i) { ret[i] *= fact_inv[i]; } return ret; } }; #line 4 "math/linalg/matrix.hpp" #include <cmath> #include <initializer_list> #include <type_traits> #line 8 "math/linalg/matrix.hpp" template <typename T> class Matrix { public: static Matrix concat(const Matrix& A, const Matrix& B) { assert(A.row == B.row); Matrix C(A.row, A.col + B.col); for (int i = 0; i < A.row; ++i) { std::ranges::copy(A[i], C[i].begin()); std::ranges::copy(B[i], C[i].begin() + A.col); } return C; } static Matrix I(int n) { Matrix ret(n); for (int i = 0; i < n; ++i) ret[i][i] = 1; return ret; } Matrix() = default; Matrix(int n) : Matrix(n, n) {} Matrix(int row, int col) : row(row), col(col), mat(row, std::vector<T>(col)) {} Matrix(const std::vector<std::vector<T>>& mat) : row(mat.size()), col(mat[0].size()), mat(mat) {} int row_size() const { return row; } int col_size() const { return col; } std::pair<int, int> shape() const { return {row, col}; } const std::vector<T>& operator[](int i) const { return mat[i]; } std::vector<T>& operator[](int i) { return mat[i]; } Matrix submatrix(int i0, int i1, int j0, int j1) const { Matrix ret(i1 - i0, j1 - j0); for (int i = i0; i < i1; ++i) { std::ranges::copy(mat[i].begin() + j0, mat[i].begin() + j1, ret.mat[i - i0].begin()); } return ret; } // --- binary operations with matrix --- Matrix& operator+=(const Matrix& rhs) { assert(shape() == rhs.shape()); for (int i = 0; i < row; ++i) { for (int j = 0; j < col; ++j) { mat[i][j] += rhs[i][j]; } } return *this; } Matrix& operator-=(const Matrix& rhs) { assert(shape() == rhs.shape()); for (int i = 0; i < row; ++i) { for (int j = 0; j < col; ++j) { mat[i][j] -= rhs[i][j]; } } return *this; } Matrix& operator*=(const Matrix& rhs) { assert(col == rhs.row); Matrix res(row, rhs.col); for (int i = 0; i < row; ++i) { for (int k = 0; k < col; ++k) { for (int j = 0; j < rhs.col; ++j) { res[i][j] += mat[i][k] * rhs.mat[k][j]; } } } return *this = res; } Matrix operator+(const Matrix& rhs) const { return Matrix(*this) += rhs; } Matrix operator-(const Matrix& rhs) const { return Matrix(*this) -= rhs; } Matrix operator*(const Matrix& rhs) const { return Matrix(*this) *= rhs; } constexpr bool operator==(const Matrix& rhs) const { if (shape() != rhs.shape()) return false; for (int i = 0; i < row; ++i) { for (int j = 0; j < col; ++j) { if (!eq(mat[i][j], rhs.mat[i][j])) return false; } } return true; } // --- scalar multiplication --- Matrix& operator*=(const T& rhs) { for (auto& row : mat) { for (auto& x : row) x *= rhs; } return *this; } Matrix& operator/=(const T& rhs) { return *this /= rhs; } Matrix operator*(const T& rhs) const { return Matrix(*this) *= rhs; } Matrix operator/(const T& rhs) const { return Matrix(*this) /= rhs; } // --- other operations for general matrices --- Matrix operator-() const { Matrix ret(*this); for (auto& row : ret.mat) { for (auto& x : row) x = -x; } return ret; } Matrix transpose() const { Matrix ret(col, row); for (int i = 0; i < col; ++i) { for (int j = 0; j < row; ++j) { ret[i][j] = mat[j][i]; } } return ret; } void reduce() { int pivot = 0; for (int j = 0; j < col; ++j) { int i = pivot; while (i < row && eq(mat[i][j], T(0))) ++i; if (i == row) continue; if (i != pivot) mat[i].swap(mat[pivot]); T pinv = T(1) / mat[pivot][j]; for (int l = j; l < col; ++l) mat[pivot][l] *= pinv; for (int k = 0; k < row; ++k) { if (k == pivot) continue; T v = mat[k][j]; for (int l = j; l < col; ++l) { mat[k][l] -= mat[pivot][l] * v; } } ++pivot; } } int rank() const { auto A = *this; A.reduce(); for (int i = 0; i < row; ++i) { bool nonzero = false; for (int j = 0; j < col; ++j) { if (!eq(A[i][j], T(0))) { nonzero = true; break; } } if (!nonzero) return i; } return row; } // --- other operations for square matrices --- void assert_square() const { assert(row == col); } Matrix pow(long long k) const { assert_square(); auto ret = I(row); auto A = *this; while (k > 0) { if (k & 1) ret *= A; A *= A; k >>= 1; } return ret; } T det() const { assert_square(); auto A = *this; T ret = 1; for (int j = 0; j < col; ++j) { int i = j; while (i < row && eq(A[i][j], T(0))) ++i; if (i == row) return 0; if (i != j) { A[i].swap(A[j]); ret = -ret; } T p = A[j][j]; ret *= p; auto pinv = T(1) / p; for (int l = j; l < col; ++l) A[j][l] *= pinv; for (int k = j + 1; k < row; ++k) { T v = A[k][j]; for (int l = j; l < col; ++l) { A[k][l] -= A[j][l] * v; } } } return ret; } Matrix inv() const { assert_square(); auto IB = concat(*this, I(row)); IB.reduce(); assert(IB.submatrix(0, row, 0, col) == I(row)); return IB.submatrix(0, row, col, 2 * col); } template <typename U, typename std::enable_if<std::is_floating_point<U>::value>::type* = nullptr> static constexpr bool eq(U a, U b) { return std::abs(a - b) < 1e-8; } template <typename U, typename std::enable_if<!std::is_floating_point< U>::value>::type* = nullptr> static constexpr bool eq(U a, U b) { return a == b; } protected: int row, col; std::vector<std::vector<T>> mat; }; #line 7 "math/linalg/characteristic_polynomial.hpp" template <typename mint> Polynomial<mint> characteristic_polynomial(Matrix<mint> mat) { mat.assert_square(); const int n = mat.row_size(); if (n == 0) return {1}; // stage 1: reduce mat to upper Hessenberg form for (int j = 0; j < n; ++j) { int i = j + 1; while (i < n && mat[i][j] == 0) ++i; if (i == n) continue; if (i != j + 1) { // swap mat[i], mat[j+1] mat[i].swap(mat[j + 1]); // swap mat[:,i], mat[:,j+1] for (int k = 0; k < n; ++k) { std::swap(mat[k][i], mat[k][j + 1]); } } mint inv = mint(1) / mat[j + 1][j]; for (int k = j + 2; k < n; ++k) { mint v = mat[k][j] * inv; // mat[k] -= mat[j+1] * v for (int l = j; l < n; ++l) { mat[k][l] -= mat[j + 1][l] * v; } // mat[:,j+1] += mat[:,k] * v for (int l = 0; l < n; ++l) { mat[l][j + 1] += mat[l][k] * v; } } } // stage 2: recursively compute char polys of leading principal submatrices std::vector<Polynomial<mint>> p(n + 1); p[0] = {1}; for (int i = 0; i < n; ++i) { p[i + 1].resize(i + 1 + 1); // (x - mat[i,i]) * p[i] for (int j = 0; j <= i; ++j) { p[i + 1][j] -= mat[i][i] * p[i][j]; p[i + 1][j + 1] += p[i][j]; } mint beta = 1; for (int k = i - 1; k >= 0; --k) { beta *= mat[k + 1][k]; p[i + 1] -= p[k] * mat[k][i] * beta; } } return p[n]; }