verify/yosupo-characteristic-polynomial.test.cpp

• View this file on GitHub
• Last update: 2026-01-05 14:20:42+00:00
• Problem: https://judge.yosupo.jp/problem/characteristic_polynomial

Depends on

• determinant-of-linear-matrix-polynomial (math/matrix/determinant-of-linear-matrix-polynomial.hpp)

Code

// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/characteristic_polynomial

#include <cstdint>
#include <iostream>
#include <vector>

#include "../math/matrix/determinant-of-linear-matrix-polynomial.hpp"

class modint998244353 {
  public:
    static constexpr std::uint32_t mod = 998244353;

    modint998244353() : v_(0) {}
    modint998244353(long long x) {
        long long y = x % static_cast<long long>(mod);
        if (y < 0)
            y += mod;
        v_ = static_cast<std::uint32_t>(y);
    }

    std::uint32_t val() const { return v_; }

    modint998244353 inv() const { return pow(*this, mod - 2); }

    modint998244353 &operator+=(const modint998244353 &rhs) {
        std::uint32_t x = v_ + rhs.v_;
        if (x >= mod)
            x -= mod;
        v_ = x;
        return *this;
    }
    modint998244353 &operator-=(const modint998244353 &rhs) {
        std::uint32_t x = (v_ >= rhs.v_) ? (v_ - rhs.v_) : (v_ + mod - rhs.v_);
        v_ = x;
        return *this;
    }
    modint998244353 &operator*=(const modint998244353 &rhs) {
        std::uint64_t x = static_cast<std::uint64_t>(v_) * rhs.v_ % mod;
        v_ = static_cast<std::uint32_t>(x);
        return *this;
    }
    modint998244353 &operator/=(const modint998244353 &rhs) {
        return *this *= rhs.inv();
    }

    friend modint998244353 operator+(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs += rhs;
    }
    friend modint998244353 operator-(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs -= rhs;
    }
    friend modint998244353 operator*(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs *= rhs;
    }
    friend modint998244353 operator/(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs /= rhs;
    }

    friend bool operator==(const modint998244353 &lhs,
                           const modint998244353 &rhs) {
        return lhs.v_ == rhs.v_;
    }
    friend bool operator!=(const modint998244353 &lhs,
                           const modint998244353 &rhs) {
        return lhs.v_ != rhs.v_;
    }

  private:
    static modint998244353 pow(modint998244353 a, long long e) {
        modint998244353 r(1);
        while (e > 0) {
            if (e & 1)
                r *= a;
            a *= a;
            e >>= 1;
        }
        return r;
    }

    std::uint32_t v_;
};

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int N;
    std::cin >> N;

    std::vector<std::vector<modint998244353>> A(
        N, std::vector<modint998244353>(N));
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            int x;
            std::cin >> x;
            A[i][j] = modint998244353(x);
        }
    }

    const std::vector<modint998244353> p = characteristic_polynomial(A);

    for (int i = 0; i < static_cast<int>(p.size()); ++i) {
        if (i)
            std::cout << ' ';
        std::cout << p[i].val();
    }
    std::cout << '\n';
    return 0;
}

#line 1 "verify/yosupo-characteristic-polynomial.test.cpp"
// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/characteristic_polynomial

#include <cstdint>
#include <iostream>
#include <vector>

#line 1 "math/matrix/determinant-of-linear-matrix-polynomial.hpp"



// 一次行列多項式 det(M0 + x M1) を係数列として求める。
// 体上でのみ動作する（除算が必要）。M0, M1 は N×N。
// 多項式は a[0] + a[1]x + ... + a[N]x^N の昇順で返す。
// M1 を掃き出して I にし、det(xI + A) を特性多項式に帰着する。
// M1 が特異でも列に x を掛ける操作を挟むことで次数 1 を保つ。
// 計算量 O(N^3)。

#include <cassert>
#include <utility>
#line 14 "math/matrix/determinant-of-linear-matrix-polynomial.hpp"

template <class T>
void hessenberg_reduction(std::vector<std::vector<T>> &matrix) {
    const int n = static_cast<int>(matrix.size());
    assert(n == 0 || static_cast<int>(matrix[0].size()) == n);
    for (int i = 1; i < n; ++i)
        assert(static_cast<int>(matrix[i].size()) == n);

    for (int r = 0; r < n - 2; ++r) {
        int piv = -1;
        for (int h = r + 1; h < n; ++h) {
            if (matrix[h][r] != T()) {
                piv = h;
                break;
            }
        }
        if (piv < 0)
            continue;

        if (piv != r + 1) {
            matrix[r + 1].swap(matrix[piv]);
            for (int i = 0; i < n; ++i)
                std::swap(matrix[i][r + 1], matrix[i][piv]);
        }

        const T rinv = T(1) / matrix[r + 1][r];
        for (int i = r + 2; i < n; ++i) {
            const T coef = matrix[i][r] * rinv;
            if (coef == T())
                continue;
            for (int j = 0; j < n; ++j)
                matrix[i][j] -= matrix[r + 1][j] * coef;
            for (int j = 0; j < n; ++j)
                matrix[j][r + 1] += matrix[j][i] * coef;
        }
    }
}

template <class T>
std::vector<T> characteristic_polynomial(std::vector<std::vector<T>> matrix) {
    const int n = static_cast<int>(matrix.size());
    assert(n == 0 || static_cast<int>(matrix[0].size()) == n);
    for (int i = 1; i < n; ++i)
        assert(static_cast<int>(matrix[i].size()) == n);

    hessenberg_reduction(matrix);

    // p[i] = det(x I_i - matrix[0..i-1][0..i-1])（係数は昇順）
    std::vector<std::vector<T>> p(n + 1);
    p[0] = {T(1)};
    for (int i = 0; i < n; ++i) {
        p[i + 1].assign(i + 2, T());

        for (int j = 0; j <= i; ++j)
            p[i + 1][j + 1] += p[i][j];
        for (int j = 0; j <= i; ++j)
            p[i + 1][j] -= p[i][j] * matrix[i][i];

        T betas = T(1);
        for (int j = i - 1; j >= 0; --j) {
            betas *= matrix[j + 1][j];
            const T hb = (T() - matrix[j][i]) * betas;
            for (int k = 0; k <= j; ++k)
                p[i + 1][k] += hb * p[j][k];
        }
    }
    return p[n];
}

template <class T>
std::vector<T>
determinant_of_linear_matrix_polynomial(std::vector<std::vector<T>> M0,
                                        std::vector<std::vector<T>> M1) {
    const int n = static_cast<int>(M0.size());
    assert(static_cast<int>(M1.size()) == n);
    if (n == 0)
        return {T(1)};
    for (int i = 0; i < n; ++i) {
        assert(static_cast<int>(M0[i].size()) == n);
        assert(static_cast<int>(M1[i].size()) == n);
    }

    int multiply_by_x = 0; // 「特定の列に x を掛ける」操作の回数
    T det_inv = T(1);      // 1 / (det A det B)

    for (int p = 0; p < n; ++p) {
        int pivot = -1;
        for (int row = p; row < n; ++row) {
            if (M1[row][p] != T()) {
                pivot = row;
                break;
            }
        }

        if (pivot < 0) {
            ++multiply_by_x;
            if (multiply_by_x > n)
                return std::vector<T>(n + 1, T());

            // x^2 の項を発生させないため、M1[0..p-1][p] を先に消す（列基本変形）。
            for (int row = 0; row < p; ++row) {
                const T v = M1[row][p];
                if (v == T())
                    continue;
                M1[row][p] = T();
                for (int i = 0; i < n; ++i)
                    M0[i][p] -= v * M0[i][row];
            }

            // (M0 + x M1) の p 列に x を掛ける（M1 の p 列が 0 のとき swap で実現できる）。
            for (int i = 0; i < n; ++i)
                std::swap(M0[i][p], M1[i][p]);

            --p; // 同じ列をやり直す（高々 n 回）
            continue;
        }

        if (pivot != p) {
            M0[pivot].swap(M0[p]);
            M1[pivot].swap(M1[p]);
            det_inv = T() - det_inv; // *= -1
        }

        const T v = M1[p][p];
        assert(v != T());
        det_inv *= v;
        const T vinv = T(1) / v;
        for (int col = 0; col < n; ++col) {
            M0[p][col] *= vinv;
            M1[p][col] *= vinv;
        }

        for (int row = 0; row < n; ++row) {
            if (row == p)
                continue;
            const T coef = M1[row][p];
            if (coef == T())
                continue;
            for (int col = 0; col < n; ++col) {
                M0[row][col] -= M0[p][col] * coef;
                M1[row][col] -= M1[p][col] * coef;
            }
        }
    }

    // M1 = I なので det(x I + M0) を求める（特性多項式 det(x I - (-M0))）。
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j)
            M0[i][j] = T() - M0[i][j];
    }
    std::vector<T> poly = characteristic_polynomial(M0);
    for (T &c : poly)
        c *= det_inv;

    if (multiply_by_x > 0)
        poly.erase(poly.begin(), poly.begin() + multiply_by_x);
    poly.resize(n + 1, T());
    return poly;
}


#line 8 "verify/yosupo-characteristic-polynomial.test.cpp"

class modint998244353 {
  public:
    static constexpr std::uint32_t mod = 998244353;

    modint998244353() : v_(0) {}
    modint998244353(long long x) {
        long long y = x % static_cast<long long>(mod);
        if (y < 0)
            y += mod;
        v_ = static_cast<std::uint32_t>(y);
    }

    std::uint32_t val() const { return v_; }

    modint998244353 inv() const { return pow(*this, mod - 2); }

    modint998244353 &operator+=(const modint998244353 &rhs) {
        std::uint32_t x = v_ + rhs.v_;
        if (x >= mod)
            x -= mod;
        v_ = x;
        return *this;
    }
    modint998244353 &operator-=(const modint998244353 &rhs) {
        std::uint32_t x = (v_ >= rhs.v_) ? (v_ - rhs.v_) : (v_ + mod - rhs.v_);
        v_ = x;
        return *this;
    }
    modint998244353 &operator*=(const modint998244353 &rhs) {
        std::uint64_t x = static_cast<std::uint64_t>(v_) * rhs.v_ % mod;
        v_ = static_cast<std::uint32_t>(x);
        return *this;
    }
    modint998244353 &operator/=(const modint998244353 &rhs) {
        return *this *= rhs.inv();
    }

    friend modint998244353 operator+(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs += rhs;
    }
    friend modint998244353 operator-(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs -= rhs;
    }
    friend modint998244353 operator*(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs *= rhs;
    }
    friend modint998244353 operator/(modint998244353 lhs,
                                     const modint998244353 &rhs) {
        return lhs /= rhs;
    }

    friend bool operator==(const modint998244353 &lhs,
                           const modint998244353 &rhs) {
        return lhs.v_ == rhs.v_;
    }
    friend bool operator!=(const modint998244353 &lhs,
                           const modint998244353 &rhs) {
        return lhs.v_ != rhs.v_;
    }

  private:
    static modint998244353 pow(modint998244353 a, long long e) {
        modint998244353 r(1);
        while (e > 0) {
            if (e & 1)
                r *= a;
            a *= a;
            e >>= 1;
        }
        return r;
    }

    std::uint32_t v_;
};

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int N;
    std::cin >> N;

    std::vector<std::vector<modint998244353>> A(
        N, std::vector<modint998244353>(N));
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            int x;
            std::cin >> x;
            A[i][j] = modint998244353(x);
        }
    }

    const std::vector<modint998244353> p = characteristic_polynomial(A);

    for (int i = 0; i < static_cast<int>(p.size()); ++i) {
        if (i)
            std::cout << ' ';
        std::cout << p[i].val();
    }
    std::cout << '\n';
    return 0;
}

Test cases

Env	Name	Status	Elapsed	Memory
g++	example_00	AC	5 ms	4 MB
g++	example_01	AC	4 ms	4 MB
g++	example_02	AC	4 ms	4 MB
g++	max_det_zero_random_00	AC	3040 ms	6 MB
g++	max_det_zero_random_01	AC	3046 ms	6 MB
g++	max_det_zero_random_02	AC	3037 ms	6 MB
g++	max_det_zero_random_03	AC	3041 ms	6 MB
g++	max_random_00	AC	3052 ms	6 MB
g++	max_random_01	AC	3061 ms	6 MB
g++	max_random_02	AC	3050 ms	6 MB
g++	max_random_03	AC	3051 ms	6 MB
g++	max_zero_00	AC	351 ms	6 MB
g++	nontrivial_frobenius_form_00	AC	346 ms	6 MB
g++	nontrivial_frobenius_form_01	AC	1500 ms	6 MB
g++	nontrivial_frobenius_form_02	AC	2271 ms	6 MB
g++	nontrivial_frobenius_form_03	AC	2577 ms	6 MB
g++	nontrivial_frobenius_form_04	AC	2742 ms	6 MB
g++	nontrivial_frobenius_form_05	AC	3045 ms	6 MB
g++	nontrivial_frobenius_form_06	AC	2987 ms	6 MB
g++	nontrivial_frobenius_form_07	AC	2914 ms	6 MB
g++	nontrivial_frobenius_form_08	AC	2891 ms	6 MB
g++	nontrivial_frobenius_form_09	AC	2849 ms	6 MB
g++	small_multiple_root_00	AC	4 ms	4 MB
g++	small_multiple_root_01	AC	4 ms	4 MB
g++	small_random_00	AC	6 ms	4 MB
g++	small_random_01	AC	4 ms	4 MB
g++	small_random_02	AC	4 ms	4 MB
g++	small_random_03	AC	4 ms	4 MB
clang++	example_00	AC	4 ms	4 MB
clang++	example_01	AC	4 ms	4 MB
clang++	example_02	AC	4 ms	4 MB
clang++	max_det_zero_random_00	AC	3438 ms	6 MB
clang++	max_det_zero_random_01	AC	3431 ms	6 MB
clang++	max_det_zero_random_02	AC	3425 ms	6 MB
clang++	max_det_zero_random_03	AC	3437 ms	6 MB
clang++	max_random_00	AC	3440 ms	6 MB
clang++	max_random_01	AC	3437 ms	6 MB
clang++	max_random_02	AC	3445 ms	6 MB
clang++	max_random_03	AC	3439 ms	6 MB
clang++	max_zero_00	AC	323 ms	6 MB
clang++	nontrivial_frobenius_form_00	AC	322 ms	6 MB
clang++	nontrivial_frobenius_form_01	AC	1590 ms	6 MB
clang++	nontrivial_frobenius_form_02	AC	2437 ms	6 MB
clang++	nontrivial_frobenius_form_03	AC	2776 ms	6 MB
clang++	nontrivial_frobenius_form_04	AC	2954 ms	6 MB
clang++	nontrivial_frobenius_form_05	AC	3433 ms	6 MB
clang++	nontrivial_frobenius_form_06	AC	3306 ms	6 MB
clang++	nontrivial_frobenius_form_07	AC	3188 ms	6 MB
clang++	nontrivial_frobenius_form_08	AC	3149 ms	6 MB
clang++	nontrivial_frobenius_form_09	AC	3087 ms	6 MB
clang++	small_multiple_root_00	AC	4 ms	4 MB
clang++	small_multiple_root_01	AC	4 ms	4 MB
clang++	small_random_00	AC	6 ms	4 MB
clang++	small_random_01	AC	4 ms	4 MB
clang++	small_random_02	AC	4 ms	4 MB
clang++	small_random_03	AC	4 ms	4 MB

Back to top page