#include #include #include #include #include #include #include #include #include template class ModInt { static_assert((Mod >> 31) == 0, "`Mod` must less than 2^(31)"); template static std::enable_if_t, unsigned> safe_mod(Int v) { using D = std::common_type_t; return (v %= (int)Mod) < 0 ? (D)(v + (int)Mod) : (D)v; } struct PrivateConstructor {}; static inline PrivateConstructor private_constructor{}; ModInt(PrivateConstructor, unsigned v) : v_(v) {} unsigned v_; public: static unsigned mod() { return Mod; } static ModInt from_raw(unsigned v) { return ModInt(private_constructor, v); } ModInt() : v_() {} template , int> = 0> ModInt(Int v) : v_(safe_mod(v)) {} template , int> = 0> ModInt(Int v) : v_(v % Mod) {} unsigned val() const { return v_; } ModInt operator-() const { return from_raw(v_ == 0 ? v_ : Mod - v_); } ModInt pow(long long e) const { if (e < 0) return inv().pow(-e); for (ModInt x(*this), res(from_raw(1));; x *= x) { if (e & 1) res *= x; if ((e >>= 1) == 0) return res; } } ModInt inv() const { int x1 = 1, x3 = 0, a = val(), b = Mod; while (b) { int q = a / b, x1_old = x1, a_old = a; x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q; } return from_raw(x1 < 0 ? x1 + (int)Mod : x1); } template std::enable_if_t div_by_2() const { if (v_ & 1) return from_raw((v_ + Mod) >> 1); return from_raw(v_ >> 1); } ModInt &operator+=(const ModInt &a) { if ((v_ += a.v_) >= Mod) v_ -= Mod; return *this; } ModInt &operator-=(const ModInt &a) { if ((v_ += Mod - a.v_) >= Mod) v_ -= Mod; return *this; } ModInt &operator*=(const ModInt &a) { v_ = (unsigned long long)v_ * a.v_ % Mod; return *this; } ModInt &operator/=(const ModInt &a) { return *this *= a.inv(); } friend ModInt operator+(const ModInt &a, const ModInt &b) { return ModInt(a) += b; } friend ModInt operator-(const ModInt &a, const ModInt &b) { return ModInt(a) -= b; } friend ModInt operator*(const ModInt &a, const ModInt &b) { return ModInt(a) *= b; } friend ModInt operator/(const ModInt &a, const ModInt &b) { return ModInt(a) /= b; } friend bool operator==(const ModInt &a, const ModInt &b) { return a.v_ == b.v_; } friend bool operator!=(const ModInt &a, const ModInt &b) { return a.v_ != b.v_; } friend std::istream &operator>>(std::istream &a, ModInt &b) { int v; a >> v; b.v_ = safe_mod(v); return a; } friend std::ostream &operator<<(std::ostream &a, const ModInt &b) { return a << b.val(); } }; template class Poly : public std::vector { public: using Base = std::vector; using Base::Base; int deg() const { for (int i = (int)this->size() - 1; i >= 0; --i) if ((*this)[i] != 0) return i; return -1; } Poly rev() const { const int d = deg(); Poly res(d + 1); for (int i = d; i >= 0; --i) res[i] = (*this)[d - i]; return res; } Tp lc() const { const int d = deg(); return d < 0 ? Tp() : (*this)[d]; } Poly &shrink() { this->resize(deg() + 1); return *this; } Poly monic() const { const int d = deg(); assert(d >= 0); Poly res(d + 1); const auto inv = lc().inv(); for (int i = 0; i <= d; ++i) res[i] = (*this)[i] * inv; return res; } Poly operator-() const { const int d = deg(); Poly res(d + 1); for (int i = 0; i <= d; ++i) res[i] = -(*this)[i]; res.shrink(); return res; } // O(deg(Q)deg(R)) std::array divmod(const Poly &R) const { const int degL = deg(), degR = R.deg(), degQ = degL - degR; assert(degR >= 0); if (degQ < 0) return {Poly(), *this}; Poly quo(degQ + 1), rem(*this); if (degQ >= 0) { const auto inv = R.lc().inv(); for (int i = degQ, n = degL; i >= 0; --i) if ((quo[i] = rem[n--] * inv) != 0) for (int j = 0; j <= degR; ++j) rem[i + j] -= quo[i] * R[j]; } rem.shrink(); return {quo, rem}; } Poly &operator+=(const Poly &R) { if (this->size() < R.size()) this->resize(R.size()); for (int i = 0; i < (int)R.size(); ++i) (*this)[i] += R[i]; return shrink(); } Poly &operator-=(const Poly &R) { if (this->size() < R.size()) this->resize(R.size()); for (int i = 0; i < (int)R.size(); ++i) (*this)[i] -= R[i]; return shrink(); } Poly &operator*=(const Poly &R) { const int degL = deg(), degR = R.deg(); if (degL < 0 || degR < 0) { this->clear(); return *this; } Poly res(degL + degR + 1); for (int i = 0; i <= degL; ++i) for (int j = 0; j <= degR; ++j) res[i + j] += (*this)[i] * R[j]; this->swap(res); return *this; } // O(min(deg(Q)^2,deg(Q)deg(R))) Poly &operator/=(const Poly &R) { const int degL = deg(), degR = R.deg(), degQ = degL - degR; assert(degR >= 0); if (degQ < 0) { this->clear(); return *this; } Poly quo(degQ + 1); const auto inv = R.lc().inv(); for (int i = 0; i <= degQ; ++i) { for (int j = 1; j <= std::min(i, degR); ++j) quo[degQ - i] += R[degR - j] * quo[degQ - i + j]; quo[degQ - i] = ((*this)[degL - i] - quo[degQ - i]) * inv; } this->swap(quo); return *this; } Poly &operator%=(const Poly &R) { const int degL = deg(), degR = R.deg(), degQ = degL - degR; assert(degR >= 0); const auto inv = R.lc().inv(); for (int i = degQ, n = degL; i >= 0; --i) if (const Tp res = (*this)[n--] * inv; res != 0) for (int j = 0; j <= degR; ++j) (*this)[i + j] -= res * R[j]; return shrink(); } Poly &operator<<=(int D) { if (D > 0) { this->insert(this->begin(), D, Tp()); } else if (D < 0) { if (-D < (int)this->size()) { this->erase(this->begin(), this->begin() + (-D)); } else { this->clear(); } } return shrink(); } Poly &operator>>=(int D) { return operator<<=(-D); } friend Poly operator+(const Poly &L, const Poly &R) { return Poly(L) += R; } friend Poly operator-(const Poly &L, const Poly &R) { return Poly(L) -= R; } friend Poly operator*(const Poly &L, const Poly &R) { return Poly(L) *= R; } friend Poly operator/(const Poly &L, const Poly &R) { return Poly(L) /= R; } friend Poly operator%(const Poly &L, const Poly &R) { return Poly(L) %= R; } friend Poly operator<<(const Poly &L, int D) { return Poly(L) <<= D; } friend Poly operator>>(const Poly &L, int D) { return Poly(L) >>= D; } friend std::ostream &operator<<(std::ostream &L, const Poly &R) { L << '['; const int d = R.deg(); if (d < 0) { L << '0'; } else { for (int i = 0; i <= d; ++i) { L << R[i]; if (i == 1) L << "*x"; if (i > 1) L << "*x^" << i; if (i != d) L << " + "; } } return L << ']'; } }; template class Matrix : public std::vector> { public: using Base = std::vector>; using Base::Base; int width() const { return this->empty() ? 0 : (int)(*this)[0].size(); } int height() const { return this->size(); } bool is_square() const { return width() == height(); } Matrix transpose() const { const int w = width(); const int h = height(); Matrix res(w, std::vector(h)); for (int i = 0; i < h; ++i) for (int j = 0; j < w; ++j) res[j][i] = (*this)[i][j]; return res; } std::vector apply(const std::vector &b) const { const int w = width(); const int h = height(); assert((int)b.size() == w); std::vector res(h); for (int i = 0; i < h; ++i) for (int j = 0; j < w; ++j) res[i] += (*this)[i][j] * b[j]; return res; } friend Matrix operator*(const Matrix &A, const Matrix &B) { const int wA = A.width(); const int hA = A.height(); assert(B.height() == wA); const int wB = B.width(); Matrix res(hA, std::vector(wB)); for (int i = 0; i < hA; ++i) for (int k = 0; k < wA; ++k) for (int j = 0; j < wB; ++j) res[i][j] += A[i][k] * B[k][j]; return res; } }; template class Basis { public: const int Dim; Matrix Vectors; // v_0, v_1, ... Matrix Augmented; Matrix Reduced; // upper triangular matrix, diagonal of Reduced = (1,...,1) // Augmented * Vectors = Reduced explicit Basis(int dim) : Dim(dim), Augmented(dim), Reduced(dim) {} int size() const { return Vectors.size(); } int dim() const { return Dim; } // if V is linear combination of v_0, ..., v_(k-1) then // returns coefficients (a_0, ..., a_(k-1)) s.t. -(a_0v_0 + ... + a_(k-1)v_(k-1)) = V std::optional> insert(const std::vector &V) { std::vector Aug(dim()), RV = V; for (int i = 0; i < dim(); ++i) { if (RV[i] == 0) continue; if (Reduced[i].empty()) { Aug[size()] = 1; const auto inv = RV[i].inv(); for (int j = i; j < dim(); ++j) RV[j] *= inv; for (int j = 0; j < dim(); ++j) Aug[j] *= inv; Augmented[i] = Aug; Reduced[i] = RV; Vectors.push_back(V); return {}; } const auto v = RV[i]; for (int j = i; j < dim(); ++j) RV[j] -= v * Reduced[i][j]; for (int j = 0; j < dim(); ++j) Aug[j] -= v * Augmented[i][j]; } return Aug; } // returns A s.t. A^(-1)MA is the linear transform respect to the basis Matrix transition_matrix() const { assert(size() == dim()); return Vectors.transpose(); } // returns A^(-1) s.t. A^(-1)MA is the linear transform respect to the basis Matrix inv_transition_matrix() const { assert(size() == dim()); auto res = Augmented; for (int i = dim() - 1; i > 0; --i) for (int j = i - 1; j >= 0; --j) for (int k = 0; k < dim(); ++k) res[j][k] -= Reduced[j][i] * res[i][k]; return res.transpose(); } }; // Compute the Frobenius form (rational canonical form) of a square matrix, // but the result is not always correct. template class Frobenius { public: // F_A = T^(-1)AT = diag(C_(p_0),...,C_(p_(k-1))) // where C_(p_j) is the companion matrix of monic polynomial P[j] // * minimal polynomial of A = p_0 // * characteristic polynomial of A = prod[0<=j InvT; std::vector> P; Matrix T; mutable std::optional> CharPoly; // see: // [1]: Elegia. A (Somehow) Simple (Randomized) Algorithm for Frobenius Form of a Matrix. // https://codeforces.com/blog/entry/124815 // [2]: Arne Storjohann. Algorithms for Matrix Canonical Forms. // https://cs.uwaterloo.ca/~astorjoh/diss2up.pdf explicit Frobenius(const Matrix &A) : N(A.height()) { assert(A.is_square()); for (;;) if (unsafe_try(A)) break; } bool unsafe_try(const Matrix &A) { static std::mt19937 rng(std::random_device{}()); std::uniform_int_distribution dis(0, Tp::mod() - 1); Basis B(N); Matrix A_B(N, std::vector(N)); // linear transform respect to basis B std::vector> V; // vectors for new basis P.clear(); while (B.size() < N) { std::vector R(N); for (int i = 0; i < N; ++i) R[i] = dis(rng); for (int deg = 0;; ++deg) { if (const auto c = B.insert(R)) { if (deg == 0) break; if (!P.empty() && deg > P.back().deg()) return false; P.emplace_back(c->begin() + (B.size() - deg), c->begin() + B.size()) .emplace_back(1); // One could check the remainder here V.emplace_back(Poly(c->begin(), c->begin() + (B.size() - deg)) / P.back()) .resize(N); V.back().at(B.size() - deg) = 1; for (int i = B.size() - deg; i < B.size() - 1; ++i) A_B[i + 1][i] = 1; for (int i = 0; i < B.size(); ++i) A_B[i][B.size() - 1] = -c->at(i); break; } R = A.apply(R); } } InvT = B.inv_transition_matrix(); T = B.transition_matrix(); if (P.size() <= 1) return true; Matrix C(N, std::vector(N)); for (int i = 0, j = 0; i < (int)V.size(); ++i) { C[j++] = V[i]; for (int k = P[i].deg(); --k; ++j) for (int l = 0; l < j; ++l) for (int m = 0; m < j; ++m) C[j][l] += A_B[l][m] * C[j - 1][m]; } for (int i = N - 1; i > 0; --i) for (int j = i - 1; j >= 0; --j) if (C[i][j] != 0) for (int k = 0; k < N; ++k) T[k][i] += C[i][j] * T[k][j], InvT[j][k] -= C[i][j] * InvT[i][k]; return true; } Matrix transition_matrix() const { return T; } Matrix inv_transition_matrix() const { return InvT; } Matrix frobenius_form() const { Matrix res(N, std::vector(N)); for (int i = 0, s = 0; i < (int)P.size(); s += P[i++].deg()) { for (int j = s; j < s + P[i].deg() - 1; ++j) res[j + 1][j] = 1; for (int j = s; j < s + P[i].deg(); ++j) res[j][s + P[i].deg() - 1] = -P[i][j - s]; } return res; } // returns (F_A)^e Matrix pow(long long e) const { assert(e >= 0); // returns x^e mod p auto pow_mod = [](auto &&pow_mod, long long e, const Poly &p) { if (p.deg() == 0) return Poly{}; if (e < p.deg()) return Poly{Tp(1)} << e; const auto half = pow_mod(pow_mod, e / 2, p); return ((half * half) << (e & 1)) % p; }; Matrix res(N, std::vector(N)); for (int i = 0, s = 0; i < (int)P.size(); s += P[i++].deg()) { auto c = pow_mod(pow_mod, e, P[i]); for (int j = 0; j < P[i].deg(); ++j) { for (int k = 0; k <= c.deg(); ++k) res[k + s][s + j] = c[k]; if (j + 1 < P[i].deg()) c = (c << 1) % P[i]; } } return res; } Poly charpoly() const { if (!CharPoly.has_value()) CharPoly.emplace(std::accumulate(P.begin(), P.end(), Poly{Tp(1)}, std::multiplies<>())); return CharPoly.value(); } // returns F(F_A) Matrix eval(Poly F) const { F %= charpoly(); Matrix res(N, std::vector(N)); if (F.deg() < 0) return res; for (int i = 0, s = 0; i < (int)P.size(); s += P[i++].deg()) { std::vector> pow_table(F.deg() + P[i].deg() + 1); pow_table[0] = Poly{Tp(1)}; for (int j = 1; j <= F.deg() + P[i].deg(); ++j) pow_table[j] = (pow_table[j - 1] << 1) % P[i]; std::vector> row(P[i].deg()); for (int j = 0; j <= F.deg(); ++j) for (int k = 0; k < P[i].deg(); ++k) row[k] += Poly{F[j]} * pow_table[j + k]; for (int j = 0; j < P[i].deg(); ++j) for (int k = 0; k <= row[j].deg(); ++k) res[k + s][s + j] = row[j][k]; } return res; } }; int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); using mint = ModInt<1000000007>; Matrix A(2, std::vector(2)); Matrix B(2, std::vector(2)); std::cin >> A[0][0] >> A[0][1] >> A[1][0] >> A[1][1]; std::cin >> B[0][0] >> B[0][1] >> B[1][0] >> B[1][1]; if (Frobenius(A).frobenius_form() == Frobenius(B).frobenius_form()) { std::cout << "Yes\n"; } else { std::cout << "No\n"; } return 0; }