Hashiryo's Library

This documentation is automatically generated by competitive-verifier/competitive-verifier

View the Project on GitHub hashiryo/Library

:heavy_check_mark: 行列の最小多項式 (src/LinearAlgebra/MinimalPolynomial.hpp)

MinimalPolynomial クラス

正方行列 $M$ と ベクトル $\boldsymbol{b}$ を与えて 最小多項式 $p(x)=p_0+p_1x+\cdots+p_{d-1}x^{d-1}+x^d$ を求める.
ただしここでの最小多項式とは、モニックな多項式であって \[ p(M)\boldsymbol{b} = p_0\boldsymbol{b}+p_1M\boldsymbol{b}+\cdots+p_{d-1}M^{d-1}\boldsymbol{b} + M^d\boldsymbol{b} = \boldsymbol{0} \] を満たす最小次数のものを指す.

ModInt クラスを前提にしている.

メンバ関数

関数 概要 計算量
MinimalPolynomial(M,b) コンストラクタ. 表現行列が $n\times n$ 正方行列な線型写像 $M$ ( Vector $\rightarrow$ Vector の関数 or Matrix クラス) とベクトル $\boldsymbol{b}$ ( Vector クラス ) を与えて最小多項式を求める $O(n^2+nT(n))$
ただし $T(n)$ は $M$ とベクトルの乗算にかかる時間
degree() 最小多項式の次元を返す $O(1)$
operator[](i) 最小多項式の$x^i$の係数を返す $O(1)$
pow(k) $M^k\boldsymbol{b}$ ( Vector クラス ) を返す $O(n^2\log k)$

その他関数

関数 概要 計算量
linear_map_det<mod_t>(M,n) 表現行列が $n\times n$ 正方行列な線型写像 $M$ ( Vector $\rightarrow$ Vector の関数 or Matrix クラス) の行列式 $\det M$ の値を返す $O(n^2+nT(n))$
ただし $T(n)$ は $M$ とベクトルの乗算にかかる時間

参考

https://yukicoder.me/wiki/black_box_linear_algebra

Depends on

Verified with

Code

#pragma once
#include <cassert>
#include "src/Internal/modint_traits.hpp"
#include "src/Math/berlekamp_massey.hpp"
#include "src/LinearAlgebra/Vector.hpp"
#include "src/Misc/rng.hpp"
// c s.t. (c[d] * M^d + c[d-1] * M^(d-1)  + ... + c[1] * M + c[0]) * b = 0
template <class mod_t, class LinMap> class MinimalPolynomial {
 std::vector<mod_t> poly, rev;
 size_t dg, n;
 std::vector<Vector<mod_t>> bs;
 static inline int deg(const std::vector<mod_t> &p) {
  for (int d= p.size() - 1;; d--)
   if (d < 0 || p[d] != mod_t()) return d;
 }
 static inline std::vector<mod_t> bostan_mori_msb(const std::vector<mod_t> &q, uint64_t k) {
  int d= deg(q);
  assert(d >= 0), assert(q[0] != mod_t());
  std::vector<mod_t> ret(std::max(d, 1));
  if (k == 0) return ret.back()= mod_t(1), ret;
  std::vector<mod_t> v(d + 1);
  for (int i= 0; i <= d; i+= 2)
   for (int j= 0; j <= d; j+= 2) v[(i + j) >> 1]+= q[i] * q[j];
  for (int i= 1; i <= d; i+= 2)
   for (int j= 1; j <= d; j+= 2) v[(i + j) >> 1]-= q[i] * q[j];
  auto w= bostan_mori_msb(v, k >> 1);
  for (int i= 2 * d - 1 - (k & 1); i >= d; i-= 2)
   for (int j= 0; j <= d; j+= 2) ret[i - d]+= q[j] * w[(i - j) >> 1];
  for (int i= 2 * d - 1 - !(k & 1); i >= d; i-= 2)
   for (int j= 1; j <= d; j+= 2) ret[i - d]-= q[j] * w[(i - j) >> 1];
  return ret;
 }
 std::vector<mod_t> x_pow_mod(uint64_t k) const {
  assert(k >= n);
  std::vector<mod_t> ret(n), u= bostan_mori_msb(rev, k - n + dg);
  for (int i= dg; i--;)
   for (int j= i + 1; j--;) ret[n - 1 - i]+= u[j] * rev[i - j];
  return ret;
 }
public:
 MinimalPolynomial(const LinMap &M, Vector<mod_t> b): n(b.size()), bs(n) {
  static_assert(is_modint_v<mod_t>);
  Vector<mod_t> a(n);
  for (size_t i= n; i--;) a[i]= rng(1, mod_t::mod() - 1);
  std::vector<mod_t> v((n + 1) << 1);
  for (size_t i= v.size(), j= 0;; b= M(b)) {
   if (j < n) bs[j]= b;
   if (v[j++]= dot(a, b); !(--i)) break;
  }
  rev= berlekamp_massey(v);
  for (auto &x: rev) x= -x;
  rev.insert(rev.begin(), 1), poly.assign(rev.rbegin(), rev.rend()), rev.erase(rev.begin() + (dg= deg(rev)) + 1, rev.end());
 }
 Vector<mod_t> pow(uint64_t k) const {  // M^k * b
  if (k < n) return bs[k];
  auto r= x_pow_mod(k);
  Vector<mod_t> ret= r[0] * bs[0];
  for (int i= r.size(); --i;) ret+= r[i] * bs[i];
  return ret;
 }
 const mod_t &operator[](size_t k) const { return poly[k]; }
 const auto begin() const { return poly.begin(); }
 const auto end() const { return poly.end(); }
 size_t degree() const { return dg; }
};
template <class mod_t, class LinMap> mod_t linear_map_det(const LinMap &M, int n) {
 Vector<mod_t> b(n);
 for (int i= n; i--;) b[i]= rng(1, mod_t::mod() - 1);
 std::vector<mod_t> D(n);
 for (auto &x: D) x= rng(1, mod_t::mod() - 1);
 auto f= [&](Vector<mod_t> a) {
  for (int i= n; i--;) a[i]*= D[i];
  return M(a);
 };
 mod_t ret= MinimalPolynomial(f, b)[0], den= 1;
 if (n & 1) ret= -ret;
 for (const auto &x: D) den*= x;
 return ret / den;
}
#line 2 "src/LinearAlgebra/MinimalPolynomial.hpp"
#include <cassert>
#line 2 "src/Internal/modint_traits.hpp"
#include <type_traits>
namespace math_internal {
struct m_b {};
struct s_b: m_b {};
}
template <class mod_t> constexpr bool is_modint_v= std::is_base_of_v<math_internal::m_b, mod_t>;
template <class mod_t> constexpr bool is_staticmodint_v= std::is_base_of_v<math_internal::s_b, mod_t>;
#line 2 "src/Math/berlekamp_massey.hpp"
#include <vector>
// a[n] = c[0] * a[n-1] + c[1] * a[n-2] + ... + c[d-1] * a[n-d]
// return c
template <class K> std::vector<K> berlekamp_massey(const std::vector<K> &a) {
 size_t n= a.size(), d= 0, m= 0, i, j;
 if (n == 0) return {};
 std::vector<K> c(n), b(n), tmp;
 K x= 1, y, coef;
 for (c[0]= b[0]= 1, i= 0; i < n; ++i) {
  for (++m, y= a[i], j= 1; j <= d; ++j) y+= c[j] * a[i - j];
  if (y == K()) continue;
  for (tmp= c, coef= y / x, j= m; j < n; ++j) c[j]-= coef * b[j - m];
  if (2 * d <= i) d= i + 1 - d, b= tmp, x= y, m= 0;
 }
 c.resize(d + 1), c.erase(c.begin());
 for (auto &x: c) x= -x;
 return c;
}
#line 2 "src/LinearAlgebra/Vector.hpp"
#include <cstdint>
#include <iostream>
#include <valarray>
namespace _la_internal {
using namespace std;
template <class R> struct Vector {
 valarray<R> dat;
 Vector()= default;
 Vector(size_t n): dat(n) {}
 Vector(size_t n, const R &v): dat(v, n) {}
 Vector(const initializer_list<R> &v): dat(v) {}
 R &operator[](int i) { return dat[i]; }
 const R &operator[](int i) const { return dat[i]; }
 bool operator==(const Vector &r) const {
  if (dat.size() != r.dat.size()) return false;
  for (int i= dat.size(); i--;)
   if (dat[i] != r.dat[i]) return false;
  return true;
 }
 bool operator!=(const Vector &r) const { return !(*this == r); }
 explicit operator bool() const { return dat.size(); }
 Vector operator-() const { return Vector(dat.size())-= *this; }
 Vector &operator+=(const Vector &r) { return dat+= r.dat, *this; }
 Vector &operator-=(const Vector &r) { return dat-= r.dat, *this; }
 Vector &operator*=(const R &r) { return dat*= r, *this; }
 Vector operator+(const Vector &r) const { return Vector(*this)+= r; }
 Vector operator-(const Vector &r) const { return Vector(*this)-= r; }
 Vector operator*(const R &r) const { return Vector(*this)*= r; }
 size_t size() const { return dat.size(); }
 friend R dot(const Vector<R> &a, const Vector<R> &b) { return assert(a.size() == b.size()), (a.dat * b.dat).sum(); }
};
using u128= __uint128_t;
using u64= uint64_t;
using u8= uint8_t;
class Ref {
 u128 *ref;
 u8 i;
public:
 Ref(u128 *ref, u8 i): ref(ref), i(i) {}
 Ref &operator=(const Ref &r) { return *this= bool(r); }
 Ref &operator=(bool b) { return *ref&= ~(u128(1) << i), *ref|= u128(b) << i, *this; }
 Ref &operator|=(bool b) { return *ref|= u128(b) << i, *this; }
 Ref &operator&=(bool b) { return *ref&= ~(u128(!b) << i), *this; }
 Ref &operator^=(bool b) { return *ref^= u128(b) << i, *this; }
 operator bool() const { return (*ref >> i) & 1; }
};
template <> class Vector<bool> {
 size_t n;
public:
 valarray<u128> dat;
 Vector(): n(0) {}
 Vector(size_t n): n(n), dat((n + 127) >> 7) {}
 Vector(size_t n, bool b): n(n), dat(-u128(b), (n + 127) >> 7) {
  if (int k= n & 127; k) dat[dat.size() - 1]&= (u128(1) << k) - 1;
 }
 Vector(const initializer_list<bool> &v): n(v.size()), dat((n + 127) >> 7) {
  int i= 0;
  for (bool b: v) dat[i >> 7]|= u128(b) << (i & 127), ++i;
 }
 Ref operator[](int i) { return {begin(dat) + (i >> 7), u8(i & 127)}; }
 bool operator[](int i) const { return (dat[i >> 7] >> (i & 127)) & 1; }
 bool operator==(const Vector &r) const {
  if (dat.size() != r.dat.size()) return false;
  for (int i= dat.size(); i--;)
   if (dat[i] != r.dat[i]) return false;
  return true;
 }
 bool operator!=(const Vector &r) const { return !(*this == r); }
 explicit operator bool() const { return n; }
 Vector operator-() const { return Vector(*this); }
 Vector &operator+=(const Vector &r) { return dat^= r.dat, *this; }
 Vector &operator-=(const Vector &r) { return dat^= r.dat, *this; }
 Vector &operator*=(bool b) { return dat*= b, *this; }
 Vector operator+(const Vector &r) const { return Vector(*this)+= r; }
 Vector operator-(const Vector &r) const { return Vector(*this)-= r; }
 Vector operator*(bool b) const { return Vector(*this)*= b; }
 size_t size() const { return n; }
 friend bool dot(const Vector<bool> &a, const Vector<bool> &b) {
  assert(a.size() == b.size());
  u128 v= 0;
  for (int i= a.dat.size(); i--;) v^= a.dat[i] & b.dat[i];
  return __builtin_parityll(v >> 64) ^ __builtin_parityll(u64(v));
 }
};
template <class R> Vector<R> operator*(const R &r, const Vector<R> &v) { return v * r; }
template <class R> ostream &operator<<(ostream &os, const Vector<R> &v) {
 os << '[';
 for (int _= 0, __= v.size(); _ < __; ++_) os << (_ ? ", " : "") << v[_];
 return os << ']';
}
}
using _la_internal::Vector;
#line 2 "src/Misc/rng.hpp"
#include <random>
#line 4 "src/Misc/rng.hpp"
uint64_t rng() {
 static uint64_t x= 10150724397891781847ULL * std::random_device{}();
 return x^= x << 7, x^= x >> 9;
}
uint64_t rng(uint64_t lim) { return rng() % lim; }
int64_t rng(int64_t l, int64_t r) { return l + rng() % (r - l); }
#line 7 "src/LinearAlgebra/MinimalPolynomial.hpp"
// c s.t. (c[d] * M^d + c[d-1] * M^(d-1)  + ... + c[1] * M + c[0]) * b = 0
template <class mod_t, class LinMap> class MinimalPolynomial {
 std::vector<mod_t> poly, rev;
 size_t dg, n;
 std::vector<Vector<mod_t>> bs;
 static inline int deg(const std::vector<mod_t> &p) {
  for (int d= p.size() - 1;; d--)
   if (d < 0 || p[d] != mod_t()) return d;
 }
 static inline std::vector<mod_t> bostan_mori_msb(const std::vector<mod_t> &q, uint64_t k) {
  int d= deg(q);
  assert(d >= 0), assert(q[0] != mod_t());
  std::vector<mod_t> ret(std::max(d, 1));
  if (k == 0) return ret.back()= mod_t(1), ret;
  std::vector<mod_t> v(d + 1);
  for (int i= 0; i <= d; i+= 2)
   for (int j= 0; j <= d; j+= 2) v[(i + j) >> 1]+= q[i] * q[j];
  for (int i= 1; i <= d; i+= 2)
   for (int j= 1; j <= d; j+= 2) v[(i + j) >> 1]-= q[i] * q[j];
  auto w= bostan_mori_msb(v, k >> 1);
  for (int i= 2 * d - 1 - (k & 1); i >= d; i-= 2)
   for (int j= 0; j <= d; j+= 2) ret[i - d]+= q[j] * w[(i - j) >> 1];
  for (int i= 2 * d - 1 - !(k & 1); i >= d; i-= 2)
   for (int j= 1; j <= d; j+= 2) ret[i - d]-= q[j] * w[(i - j) >> 1];
  return ret;
 }
 std::vector<mod_t> x_pow_mod(uint64_t k) const {
  assert(k >= n);
  std::vector<mod_t> ret(n), u= bostan_mori_msb(rev, k - n + dg);
  for (int i= dg; i--;)
   for (int j= i + 1; j--;) ret[n - 1 - i]+= u[j] * rev[i - j];
  return ret;
 }
public:
 MinimalPolynomial(const LinMap &M, Vector<mod_t> b): n(b.size()), bs(n) {
  static_assert(is_modint_v<mod_t>);
  Vector<mod_t> a(n);
  for (size_t i= n; i--;) a[i]= rng(1, mod_t::mod() - 1);
  std::vector<mod_t> v((n + 1) << 1);
  for (size_t i= v.size(), j= 0;; b= M(b)) {
   if (j < n) bs[j]= b;
   if (v[j++]= dot(a, b); !(--i)) break;
  }
  rev= berlekamp_massey(v);
  for (auto &x: rev) x= -x;
  rev.insert(rev.begin(), 1), poly.assign(rev.rbegin(), rev.rend()), rev.erase(rev.begin() + (dg= deg(rev)) + 1, rev.end());
 }
 Vector<mod_t> pow(uint64_t k) const {  // M^k * b
  if (k < n) return bs[k];
  auto r= x_pow_mod(k);
  Vector<mod_t> ret= r[0] * bs[0];
  for (int i= r.size(); --i;) ret+= r[i] * bs[i];
  return ret;
 }
 const mod_t &operator[](size_t k) const { return poly[k]; }
 const auto begin() const { return poly.begin(); }
 const auto end() const { return poly.end(); }
 size_t degree() const { return dg; }
};
template <class mod_t, class LinMap> mod_t linear_map_det(const LinMap &M, int n) {
 Vector<mod_t> b(n);
 for (int i= n; i--;) b[i]= rng(1, mod_t::mod() - 1);
 std::vector<mod_t> D(n);
 for (auto &x: D) x= rng(1, mod_t::mod() - 1);
 auto f= [&](Vector<mod_t> a) {
  for (int i= n; i--;) a[i]*= D[i];
  return M(a);
 };
 mod_t ret= MinimalPolynomial(f, b)[0], den= 1;
 if (n & 1) ret= -ret;
 for (const auto &x: D) den*= x;
 return ret / den;
}
Back to top page