verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp

• View this file on GitHub
• Last update: 2026-07-11 16:54:42+09:00

Depends on

• 正方行列の 1 次式の行列式 (yukicoder No.1907 DETERMINATION) (math/matrix/determinant-of-linear-matrix-polynomial.hpp)

Code

// competitive-verifier: STANDALONE

#include <cassert>
#include <vector>

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

namespace {
struct ModInt101 {
    static constexpr int mod = 101;
    int value;

    ModInt101(long long value = 0) {
        value %= mod;
        if (value < 0) {
            value += mod;
        }
        this->value = static_cast<int>(value);
    }

    ModInt101 &operator+=(const ModInt101 &other) {
        value += other.value;
        if (value >= mod) {
            value -= mod;
        }
        return *this;
    }

    ModInt101 &operator-=(const ModInt101 &other) {
        value -= other.value;
        if (value < 0) {
            value += mod;
        }
        return *this;
    }

    ModInt101 &operator*=(const ModInt101 &other) {
        value = value * other.value % mod;
        return *this;
    }

    ModInt101 &operator/=(const ModInt101 &other) {
        assert(other != ModInt101());
        return *this *= power(other, mod - 2);
    }

    friend ModInt101 operator-(ModInt101 lhs, const ModInt101 &rhs) {
        return lhs -= rhs;
    }

    friend ModInt101 operator*(ModInt101 lhs, const ModInt101 &rhs) {
        return lhs *= rhs;
    }

    friend ModInt101 operator/(ModInt101 lhs, const ModInt101 &rhs) {
        return lhs /= rhs;
    }

    friend bool operator==(const ModInt101 &lhs, const ModInt101 &rhs) {
        return lhs.value == rhs.value;
    }

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

  private:
    static ModInt101 power(ModInt101 base, int exponent) {
        ModInt101 result(1);
        while (exponent > 0) {
            if (exponent % 2 == 1) {
                result *= base;
            }
            base *= base;
            exponent /= 2;
        }
        return result;
    }
};

using Matrix = std::vector<std::vector<ModInt101>>;

ModInt101 determinant(Matrix matrix) {
    const int n = static_cast<int>(matrix.size());
    ModInt101 result(1);
    for (int column = 0; column < n; ++column) {
        int pivot = -1;
        for (int row = column; row < n; ++row) {
            if (matrix[row][column] != ModInt101()) {
                pivot = row;
                break;
            }
        }
        if (pivot < 0) {
            return ModInt101();
        }
        if (pivot != column) {
            matrix[pivot].swap(matrix[column]);
            result = ModInt101() - result;
        }

        result *= matrix[column][column];
        const ModInt101 pivot_inverse = ModInt101(1) / matrix[column][column];
        for (int row = column + 1; row < n; ++row) {
            const ModInt101 factor = matrix[row][column] * pivot_inverse;
            for (int j = column + 1; j < n; ++j) {
                matrix[row][j] -= matrix[column][j] * factor;
            }
        }
    }
    return result;
}

ModInt101 evaluate(const std::vector<ModInt101> &polynomial, ModInt101 x) {
    ModInt101 result;
    for (int i = static_cast<int>(polynomial.size()) - 1; i >= 0; --i) {
        result *= x;
        result += polynomial[i];
    }
    return result;
}

Matrix matrix_from_mask(int n, int mask) {
    Matrix matrix(n, std::vector<ModInt101>(n));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            matrix[i][j] = ModInt101((mask >> (i * n + j)) % 2);
        }
    }
    return matrix;
}

void check_characteristic_polynomial(const Matrix &matrix) {
    const int n = static_cast<int>(matrix.size());
    const std::vector<ModInt101> polynomial = characteristic_polynomial(matrix);
    assert(static_cast<int>(polynomial.size()) == n + 1);

    Matrix hessenberg = matrix;
    hessenberg_reduction(hessenberg);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j + 1 < i; ++j) {
            assert(hessenberg[i][j] == ModInt101());
        }
    }

    for (int x = 0; x <= n; ++x) {
        Matrix shifted = matrix;
        Matrix shifted_hessenberg = hessenberg;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                const ModInt101 diagonal = i == j ? ModInt101(x) : ModInt101();
                shifted[i][j] = diagonal - shifted[i][j];
                shifted_hessenberg[i][j] = diagonal - shifted_hessenberg[i][j];
            }
        }
        const ModInt101 expected = determinant(shifted);
        assert(evaluate(polynomial, ModInt101(x)) == expected);
        assert(determinant(shifted_hessenberg) == expected);
    }
}

void check_linear_matrix_polynomial(const Matrix &matrix_0,
                                    const Matrix &matrix_1) {
    const int n = static_cast<int>(matrix_0.size());
    const std::vector<ModInt101> polynomial =
        determinant_of_linear_matrix_polynomial(matrix_0, matrix_1);
    assert(static_cast<int>(polynomial.size()) == n + 1);

    for (int x = 0; x <= n; ++x) {
        Matrix evaluated = matrix_0;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                evaluated[i][j] += ModInt101(x) * matrix_1[i][j];
            }
        }
        assert(evaluate(polynomial, ModInt101(x)) == determinant(evaluated));
    }
}

void self_test() {
    check_characteristic_polynomial({});
    check_linear_matrix_polynomial({}, {});

    for (int n = 1; n <= 3; ++n) {
        const int state_count = 1 << (n * n);
        for (int mask = 0; mask < state_count; ++mask) {
            check_characteristic_polynomial(matrix_from_mask(n, mask));
        }
    }

    for (int n = 1; n <= 2; ++n) {
        const int state_count = 1 << (n * n);
        for (int mask_0 = 0; mask_0 < state_count; ++mask_0) {
            for (int mask_1 = 0; mask_1 < state_count; ++mask_1) {
                check_linear_matrix_polynomial(matrix_from_mask(n, mask_0),
                                               matrix_from_mask(n, mask_1));
            }
        }
    }

    constexpr int state_count = 1 << 9;
    for (int mask = 0; mask < state_count; ++mask) {
        check_linear_matrix_polynomial(matrix_from_mask(3, mask),
                                       matrix_from_mask(3, 0));
        check_linear_matrix_polynomial(
            matrix_from_mask(3, mask),
            matrix_from_mask(3, (mask * 137 + 91) & (state_count - 1)));
    }
}
} // namespace

int main() {
    self_test();

    return 0;
}

#line 1 "verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp"
// competitive-verifier: STANDALONE

#include <cassert>
#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)。

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

namespace determinant_of_linear_matrix_polynomial_internal {
template <class T>
bool is_square_matrix(const std::vector<std::vector<T>> &matrix) {
    const int n = static_cast<int>(matrix.size());
    for (const std::vector<T> &row : matrix) {
        if (static_cast<int>(row.size()) != n) {
            return false;
        }
    }
    return true;
}

template <class T>
void hessenberg_reduction(std::vector<std::vector<T>> &matrix) {
    const int n = static_cast<int>(matrix.size());
    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;
            }
            matrix[i][r] = T();
            for (int j = r + 1; 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());
    determinant_of_linear_matrix_polynomial_internal::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());

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

        T betas = T(1);
        for (int j = i - 1; j >= 0; --j) {
            betas *= matrix[j + 1][j];
            if (betas == T()) {
                break;
            }
            if (matrix[j][i] == T()) {
                continue;
            }
            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];
}
} // namespace determinant_of_linear_matrix_polynomial_internal

template <class T>
void hessenberg_reduction(std::vector<std::vector<T>> &matrix) {
    assert(determinant_of_linear_matrix_polynomial_internal::is_square_matrix(
        matrix));

    determinant_of_linear_matrix_polynomial_internal::hessenberg_reduction(
        matrix);
}

template <class T>
std::vector<T> characteristic_polynomial(std::vector<std::vector<T>> matrix) {
    assert(determinant_of_linear_matrix_polynomial_internal::is_square_matrix(
        matrix));

    return determinant_of_linear_matrix_polynomial_internal::
        characteristic_polynomial(std::move(matrix));
}

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);
    assert(
        determinant_of_linear_matrix_polynomial_internal::is_square_matrix(M0));
    assert(
        determinant_of_linear_matrix_polynomial_internal::is_square_matrix(M1));
    if (n == 0) {
        return {T(1)};
    }

    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];
                }
            }

            // M1 の p 列が 0 なので、M0 の p 列を M1 に移して列に x を掛ける。
            for (int i = 0; i < n; ++i) {
                M1[i][p] = std::move(M0[i][p]);
                M0[i][p] = T();
            }

            --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];
        det_inv *= v;
        const T vinv = T(1) / v;
        for (int col = 0; col < n; ++col) {
            M0[p][col] *= vinv;
        }
        M1[p][p] = T(1);
        for (int col = p + 1; col < n; ++col) {
            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][p] = T();
            for (int col = p + 1; col < n; ++col) {
                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 = determinant_of_linear_matrix_polynomial_internal::
        characteristic_polynomial(std::move(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 7 "verify/standalone-determinant-of-linear-matrix-polynomial.test.cpp"

namespace {
struct ModInt101 {
    static constexpr int mod = 101;
    int value;

    ModInt101(long long value = 0) {
        value %= mod;
        if (value < 0) {
            value += mod;
        }
        this->value = static_cast<int>(value);
    }

    ModInt101 &operator+=(const ModInt101 &other) {
        value += other.value;
        if (value >= mod) {
            value -= mod;
        }
        return *this;
    }

    ModInt101 &operator-=(const ModInt101 &other) {
        value -= other.value;
        if (value < 0) {
            value += mod;
        }
        return *this;
    }

    ModInt101 &operator*=(const ModInt101 &other) {
        value = value * other.value % mod;
        return *this;
    }

    ModInt101 &operator/=(const ModInt101 &other) {
        assert(other != ModInt101());
        return *this *= power(other, mod - 2);
    }

    friend ModInt101 operator-(ModInt101 lhs, const ModInt101 &rhs) {
        return lhs -= rhs;
    }

    friend ModInt101 operator*(ModInt101 lhs, const ModInt101 &rhs) {
        return lhs *= rhs;
    }

    friend ModInt101 operator/(ModInt101 lhs, const ModInt101 &rhs) {
        return lhs /= rhs;
    }

    friend bool operator==(const ModInt101 &lhs, const ModInt101 &rhs) {
        return lhs.value == rhs.value;
    }

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

  private:
    static ModInt101 power(ModInt101 base, int exponent) {
        ModInt101 result(1);
        while (exponent > 0) {
            if (exponent % 2 == 1) {
                result *= base;
            }
            base *= base;
            exponent /= 2;
        }
        return result;
    }
};

using Matrix = std::vector<std::vector<ModInt101>>;

ModInt101 determinant(Matrix matrix) {
    const int n = static_cast<int>(matrix.size());
    ModInt101 result(1);
    for (int column = 0; column < n; ++column) {
        int pivot = -1;
        for (int row = column; row < n; ++row) {
            if (matrix[row][column] != ModInt101()) {
                pivot = row;
                break;
            }
        }
        if (pivot < 0) {
            return ModInt101();
        }
        if (pivot != column) {
            matrix[pivot].swap(matrix[column]);
            result = ModInt101() - result;
        }

        result *= matrix[column][column];
        const ModInt101 pivot_inverse = ModInt101(1) / matrix[column][column];
        for (int row = column + 1; row < n; ++row) {
            const ModInt101 factor = matrix[row][column] * pivot_inverse;
            for (int j = column + 1; j < n; ++j) {
                matrix[row][j] -= matrix[column][j] * factor;
            }
        }
    }
    return result;
}

ModInt101 evaluate(const std::vector<ModInt101> &polynomial, ModInt101 x) {
    ModInt101 result;
    for (int i = static_cast<int>(polynomial.size()) - 1; i >= 0; --i) {
        result *= x;
        result += polynomial[i];
    }
    return result;
}

Matrix matrix_from_mask(int n, int mask) {
    Matrix matrix(n, std::vector<ModInt101>(n));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            matrix[i][j] = ModInt101((mask >> (i * n + j)) % 2);
        }
    }
    return matrix;
}

void check_characteristic_polynomial(const Matrix &matrix) {
    const int n = static_cast<int>(matrix.size());
    const std::vector<ModInt101> polynomial = characteristic_polynomial(matrix);
    assert(static_cast<int>(polynomial.size()) == n + 1);

    Matrix hessenberg = matrix;
    hessenberg_reduction(hessenberg);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j + 1 < i; ++j) {
            assert(hessenberg[i][j] == ModInt101());
        }
    }

    for (int x = 0; x <= n; ++x) {
        Matrix shifted = matrix;
        Matrix shifted_hessenberg = hessenberg;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                const ModInt101 diagonal = i == j ? ModInt101(x) : ModInt101();
                shifted[i][j] = diagonal - shifted[i][j];
                shifted_hessenberg[i][j] = diagonal - shifted_hessenberg[i][j];
            }
        }
        const ModInt101 expected = determinant(shifted);
        assert(evaluate(polynomial, ModInt101(x)) == expected);
        assert(determinant(shifted_hessenberg) == expected);
    }
}

void check_linear_matrix_polynomial(const Matrix &matrix_0,
                                    const Matrix &matrix_1) {
    const int n = static_cast<int>(matrix_0.size());
    const std::vector<ModInt101> polynomial =
        determinant_of_linear_matrix_polynomial(matrix_0, matrix_1);
    assert(static_cast<int>(polynomial.size()) == n + 1);

    for (int x = 0; x <= n; ++x) {
        Matrix evaluated = matrix_0;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                evaluated[i][j] += ModInt101(x) * matrix_1[i][j];
            }
        }
        assert(evaluate(polynomial, ModInt101(x)) == determinant(evaluated));
    }
}

void self_test() {
    check_characteristic_polynomial({});
    check_linear_matrix_polynomial({}, {});

    for (int n = 1; n <= 3; ++n) {
        const int state_count = 1 << (n * n);
        for (int mask = 0; mask < state_count; ++mask) {
            check_characteristic_polynomial(matrix_from_mask(n, mask));
        }
    }

    for (int n = 1; n <= 2; ++n) {
        const int state_count = 1 << (n * n);
        for (int mask_0 = 0; mask_0 < state_count; ++mask_0) {
            for (int mask_1 = 0; mask_1 < state_count; ++mask_1) {
                check_linear_matrix_polynomial(matrix_from_mask(n, mask_0),
                                               matrix_from_mask(n, mask_1));
            }
        }
    }

    constexpr int state_count = 1 << 9;
    for (int mask = 0; mask < state_count; ++mask) {
        check_linear_matrix_polynomial(matrix_from_mask(3, mask),
                                       matrix_from_mask(3, 0));
        check_linear_matrix_polynomial(
            matrix_from_mask(3, mask),
            matrix_from_mask(3, (mask * 137 + 91) & (state_count - 1)));
    }
}
} // namespace

int main() {
    self_test();

    return 0;
}

Back to top page