verify/yukicoder-1907-determination.test.cpp

• View this file on GitHub
• Last update: 2026-01-05 14:20:42+00:00
• Problem: https://yukicoder.me/problems/no/1907

Depends on

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

Code

// competitive-verifier: PROBLEM https://yukicoder.me/problems/no/1907

#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>> M0(
        N, std::vector<modint998244353>(N));
    std::vector<std::vector<modint998244353>> M1(
        N, std::vector<modint998244353>(N));

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

    const std::vector<modint998244353> ans =
        determinant_of_linear_matrix_polynomial(M0, M1);
    for (const auto &a : ans) {
        std::cout << a.val() << '\n';
    }
    return 0;
}

#line 1 "verify/yukicoder-1907-determination.test.cpp"
// competitive-verifier: PROBLEM https://yukicoder.me/problems/no/1907

#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/yukicoder-1907-determination.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>> M0(
        N, std::vector<modint998244353>(N));
    std::vector<std::vector<modint998244353>> M1(
        N, std::vector<modint998244353>(N));

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

    const std::vector<modint998244353> ans =
        determinant_of_linear_matrix_polynomial(M0, M1);
    for (const auto &a : ans) {
        std::cout << a.val() << '\n';
    }
    return 0;
}

Test cases

Env	Name	Status	Elapsed	Memory
g++	01_sample_01.txt	AC	5 ms	4 MB
g++	01_sample_02.txt	AC	4 ms	4 MB
g++	01_sample_03.txt	AC	4 ms	4 MB
g++	01_sample_04.txt	AC	4 ms	3 MB
g++	02_small_random_01	AC	5 ms	3 MB
g++	02_small_random_02	AC	5 ms	4 MB
g++	02_small_random_03	AC	4 ms	4 MB
g++	03_large_random_01	AC	1717 ms	5 MB
g++	03_large_random_02	AC	681 ms	5 MB
g++	03_large_random_03	AC	1181 ms	5 MB
g++	04_large_random_01	AC	3404 ms	7 MB
g++	04_large_random_02	AC	843 ms	6 MB
g++	04_large_random_03	AC	3999 ms	7 MB
g++	10_sparse_01	AC	3606 ms	7 MB
g++	10_sparse_02	AC	3227 ms	7 MB
g++	10_sparse_03	AC	598 ms	5 MB
g++	11_dense_01	AC	223 ms	4 MB
g++	11_dense_02	AC	3069 ms	7 MB
g++	11_dense_03	AC	2406 ms	6 MB
g++	12_wide_dense_01	AC	58 ms	4 MB
g++	12_wide_dense_02	AC	3536 ms	7 MB
g++	12_wide_dense_03	AC	323 ms	4 MB
g++	21_corner_01	AC	2860 ms	6 MB
g++	21_corner_02	AC	4064 ms	7 MB
g++	21_corner_03	AC	1208 ms	5 MB
g++	31_triangle_01	AC	6 ms	4 MB
g++	31_triangle_02	AC	4256 ms	7 MB
g++	31_triangle_03	AC	4775 ms	7 MB
g++	31_triangle_04	AC	4772 ms	7 MB
g++	31_triangle_05	AC	4740 ms	7 MB
g++	41_discrete_nonzero_01	AC	6 ms	4 MB
g++	41_discrete_nonzero_02	AC	4333 ms	7 MB
g++	41_discrete_nonzero_03	AC	4220 ms	7 MB
g++	41_discrete_nonzero_04	AC	4249 ms	7 MB
g++	41_discrete_nonzero_05	AC	4256 ms	7 MB
g++	51_retry_01	AC	5 ms	4 MB
g++	51_retry_02	AC	4 ms	4 MB
g++	51_retry_03	AC	5 ms	4 MB
g++	51_retry_04	AC	4243 ms	7 MB
g++	51_retry_05	AC	4253 ms	7 MB
g++	51_retry_06	AC	5681 ms	7 MB
g++	51_retry_07	AC	4243 ms	7 MB
g++	51_retry_08	AC	5684 ms	7 MB
g++	51_retry_09	AC	5687 ms	7 MB
g++	51_retry_10	AC	4674 ms	7 MB
g++	51_retry_11	AC	4939 ms	7 MB
g++	61_notrandom_killer_01	AC	4133 ms	7 MB
g++	61_notrandom_killer_02	AC	4159 ms	7 MB
g++	61_notrandom_killer_03	AC	4197 ms	7 MB
g++	61_notrandom_killer_04	AC	4228 ms	7 MB
g++	61_notrandom_killer_05	AC	4191 ms	7 MB
g++	61_notrandom_killer_06	AC	4214 ms	7 MB
g++	91_multiply_loop_01	AC	5 ms	4 MB
g++	91_multiply_loop_02	AC	4057 ms	6 MB
g++	91_multiply_loop_03	AC	4058 ms	6 MB
g++	92_multiply_loop_01	AC	5 ms	4 MB
g++	92_multiply_loop_02	AC	2619 ms	6 MB
g++	92_multiply_loop_03	AC	2628 ms	6 MB
g++	99_fixeddeg_01	AC	2217 ms	7 MB
g++	99_fixeddeg_02	AC	2763 ms	7 MB
g++	99_fixeddeg_03	AC	2812 ms	7 MB
g++	99_fixeddeg_04	AC	3014 ms	7 MB
g++	99_fixeddeg_05	AC	2799 ms	7 MB
g++	99_fixeddeg_06	AC	4247 ms	7 MB
g++	99_fixeddeg_07	AC	4 ms	4 MB
g++	99_fixeddeg_08	AC	4 ms	4 MB
g++	99_fixeddeg_09	AC	4 ms	4 MB
clang++	01_sample_01.txt	AC	5 ms	4 MB
clang++	01_sample_02.txt	AC	4 ms	4 MB
clang++	01_sample_03.txt	AC	4 ms	4 MB
clang++	01_sample_04.txt	AC	4 ms	4 MB
clang++	02_small_random_01	AC	5 ms	4 MB
clang++	02_small_random_02	AC	5 ms	4 MB
clang++	02_small_random_03	AC	4 ms	4 MB
clang++	03_large_random_01	AC	1783 ms	5 MB
clang++	03_large_random_02	AC	706 ms	5 MB
clang++	03_large_random_03	AC	1224 ms	5 MB
clang++	04_large_random_01	AC	3569 ms	7 MB
clang++	04_large_random_02	AC	877 ms	6 MB
clang++	04_large_random_03	AC	4149 ms	7 MB
clang++	10_sparse_01	AC	3743 ms	7 MB
clang++	10_sparse_02	AC	3357 ms	7 MB
clang++	10_sparse_03	AC	625 ms	5 MB
clang++	11_dense_01	AC	231 ms	4 MB
clang++	11_dense_02	AC	3224 ms	7 MB
clang++	11_dense_03	AC	2512 ms	6 MB
clang++	12_wide_dense_01	AC	60 ms	4 MB
clang++	12_wide_dense_02	AC	3676 ms	7 MB
clang++	12_wide_dense_03	AC	335 ms	4 MB
clang++	21_corner_01	AC	3014 ms	6 MB
clang++	21_corner_02	AC	4259 ms	7 MB
clang++	21_corner_03	AC	1255 ms	5 MB
clang++	31_triangle_01	AC	5 ms	4 MB
clang++	31_triangle_02	AC	4421 ms	7 MB
clang++	31_triangle_03	AC	4994 ms	7 MB
clang++	31_triangle_04	AC	4999 ms	7 MB
clang++	31_triangle_05	AC	4964 ms	7 MB
clang++	41_discrete_nonzero_01	AC	6 ms	4 MB
clang++	41_discrete_nonzero_02	AC	4489 ms	7 MB
clang++	41_discrete_nonzero_03	AC	4387 ms	7 MB
clang++	41_discrete_nonzero_04	AC	4417 ms	7 MB
clang++	41_discrete_nonzero_05	AC	4431 ms	7 MB
clang++	51_retry_01	AC	6 ms	4 MB
clang++	51_retry_02	AC	5 ms	4 MB
clang++	51_retry_03	AC	5 ms	4 MB
clang++	51_retry_04	AC	4421 ms	7 MB
clang++	51_retry_05	AC	4417 ms	7 MB
clang++	51_retry_06	AC	5982 ms	7 MB
clang++	51_retry_07	AC	4418 ms	7 MB
clang++	51_retry_08	AC	6004 ms	7 MB
clang++	51_retry_09	AC	5986 ms	7 MB
clang++	51_retry_10	AC	4885 ms	7 MB
clang++	51_retry_11	AC	5175 ms	7 MB
clang++	61_notrandom_killer_01	AC	4301 ms	7 MB
clang++	61_notrandom_killer_02	AC	4329 ms	7 MB
clang++	61_notrandom_killer_03	AC	4367 ms	7 MB
clang++	61_notrandom_killer_04	AC	4393 ms	7 MB
clang++	61_notrandom_killer_05	AC	4369 ms	7 MB
clang++	61_notrandom_killer_06	AC	4361 ms	7 MB
clang++	91_multiply_loop_01	AC	5 ms	4 MB
clang++	91_multiply_loop_02	AC	4302 ms	6 MB
clang++	91_multiply_loop_03	AC	4307 ms	6 MB
clang++	92_multiply_loop_01	AC	5 ms	4 MB
clang++	92_multiply_loop_02	AC	2736 ms	6 MB
clang++	92_multiply_loop_03	AC	2742 ms	6 MB
clang++	99_fixeddeg_01	AC	2371 ms	7 MB
clang++	99_fixeddeg_02	AC	2900 ms	7 MB
clang++	99_fixeddeg_03	AC	2918 ms	7 MB
clang++	99_fixeddeg_04	AC	3153 ms	7 MB
clang++	99_fixeddeg_05	AC	2915 ms	7 MB
clang++	99_fixeddeg_06	AC	4421 ms	7 MB
clang++	99_fixeddeg_07	AC	5 ms	4 MB
clang++	99_fixeddeg_08	AC	4 ms	4 MB
clang++	99_fixeddeg_09	AC	4 ms	4 MB

Back to top page