#1170083 (C++17) No.3231 2×2行列相似判定 〜hard〜

• 問題: No.3231 2×2行列相似判定 〜hard〜
• ユーザー: hly1204
• 提出日時: 2026-07-01 19:58:04
• 言語: C++17 (gcc 15.2.0 + boost 1.89.0)
• 結果: AC
• 実行時間: 2 ms / 2,000 ms
• コード長: 17,181 bytes
• コンパイル時間: 1,538 ms
• コンパイル使用メモリ: 155,288 KB
• 実行使用メモリ: 6,528 KB
• 最終ジャッジ日時: 2026-07-01 19:58:09
• 合計ジャッジ時間: 2,973 ms
• ジャッジサーバーID （参考情報）: judge1_0 / judge2_0

ファイルパターン	結果
sample	AC * 3
other	AC * 37

ソースコード

#include <array>
#include <cassert>
#include <functional>
#include <iostream>
#include <numeric>
#include <optional>
#include <random>
#include <type_traits>
#include <vector>

template <unsigned Mod> class ModInt {
    static_assert((Mod >> 31) == 0, "`Mod` must less than 2^(31)");
    template <typename Int>
    static std::enable_if_t<std::is_integral_v<Int>, unsigned> safe_mod(Int v) {
        using D = std::common_type_t<Int, unsigned>;
        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 <typename Int, typename std::enable_if_t<std::is_signed_v<Int>, int> = 0> ModInt(Int v)
        : v_(safe_mod(v)) {}
    template <typename Int, typename std::enable_if_t<std::is_unsigned_v<Int>, 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 <bool Odd = (Mod & 1)> std::enable_if_t<Odd, ModInt> 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 <typename Tp> class Poly : public std::vector<Tp> {
public:
    using Base = std::vector<Tp>;
    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<Poly, 2> 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 <typename Tp> class Matrix : public std::vector<std::vector<Tp>> {
public:
    using Base = std::vector<std::vector<Tp>>;
    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<Tp>(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<Tp> apply(const std::vector<Tp> &b) const {
        const int w = width();
        const int h = height();
        assert((int)b.size() == w);
        std::vector<Tp> 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<Tp>(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 <typename Tp> class Basis {
public:
    const int Dim;
    Matrix<Tp> Vectors; // v_0, v_1, ...
    Matrix<Tp> Augmented;
    Matrix<Tp> 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<std::vector<Tp>> insert(const std::vector<Tp> &V) {
        std::vector<Tp> 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<Tp> 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<Tp> 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 <typename Tp> 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<k] p_j
    const int N;
    Matrix<Tp> InvT;
    std::vector<Poly<Tp>> P;
    Matrix<Tp> T;
    mutable std::optional<Poly<Tp>> 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<Tp> &A) : N(A.height()) {
        assert(A.is_square());
        for (;;)
            if (unsafe_try(A)) break;
    }

    bool unsafe_try(const Matrix<Tp> &A) {
        static std::mt19937 rng(std::random_device{}());
        std::uniform_int_distribution<decltype(Tp::mod())> dis(0, Tp::mod() - 1);

        Basis<Tp> B(N);
        Matrix<Tp> A_B(N, std::vector<Tp>(N)); // linear transform respect to basis B
        std::vector<std::vector<Tp>> V;        // vectors for new basis
        P.clear();

        while (B.size() < N) {
            std::vector<Tp> 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<Tp>(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<Tp> C(N, std::vector<Tp>(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<Tp> transition_matrix() const { return T; }
    Matrix<Tp> inv_transition_matrix() const { return InvT; }

    Matrix<Tp> frobenius_form() const {
        Matrix<Tp> res(N, std::vector<Tp>(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<Tp> pow(long long e) const {
        assert(e >= 0);

        // returns x^e mod p
        auto pow_mod = [](auto &&pow_mod, long long e, const Poly<Tp> &p) {
            if (p.deg() == 0) return Poly<Tp>{};
            if (e < p.deg()) return Poly<Tp>{Tp(1)} << e;
            const auto half = pow_mod(pow_mod, e / 2, p);
            return ((half * half) << (e & 1)) % p;
        };

        Matrix<Tp> res(N, std::vector<Tp>(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<Tp> charpoly() const {
        if (!CharPoly.has_value())
            CharPoly.emplace(std::accumulate(P.begin(), P.end(), Poly<Tp>{Tp(1)}, std::multiplies<>()));
        return CharPoly.value();
    }

    // returns F(F_A)
    Matrix<Tp> eval(Poly<Tp> F) const {
        F %= charpoly();
        Matrix<Tp> res(N, std::vector<Tp>(N));
        if (F.deg() < 0) return res;
        for (int i = 0, s = 0; i < (int)P.size(); s += P[i++].deg()) {
            std::vector<Poly<Tp>> pow_table(F.deg() + P[i].deg() + 1);
            pow_table[0] = Poly<Tp>{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<Poly<Tp>> row(P[i].deg());
            for (int j = 0; j <= F.deg(); ++j)
                for (int k = 0; k < P[i].deg(); ++k) row[k] += Poly<Tp>{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<mint> A(2, std::vector<mint>(2));
    Matrix<mint> B(2, std::vector<mint>(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;
}