m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: Characteristic Polynomial
(math/matrix/characteristic_polynomial.hpp)

Overview

characteristic_polynomial computes the monic characteristic polynomial of a square matrix:

\[p(x) = \det(xI - A).\] 
#include "math/matrix/characteristic_polynomial.hpp"

The function and Matrix are in m1une::matrix.

Interface

template <class T>
std::vector<T> characteristic_polynomial(Matrix<T> matrix);
Function Description Complexity
characteristic_polynomial(matrix) Returns the coefficients of $\det(xI-A)$ in ascending degree order. $O(N^3)$ time and $O(N^2)$ additional memory

The input must be square. For an N by N matrix, the returned vector has N + 1 elements and satisfies

\[p(x) = p_0 + p_1x + \dots + p_Nx^N.\]

In particular, result[N] is one and result[0] is $(-1)^N\det(A)$. The characteristic polynomial of the empty 0 by 0 matrix is the constant polynomial {1}.

The matrix is passed by value and reduced internally, so the caller’s matrix is not modified unless it is explicitly moved into the function.

T must represent a field: it must support construction from zero and one, equality, addition, subtraction, multiplication, and division by a nonzero value. Fixed-prime modular integers are the main intended type. The algorithm uses exact zero comparisons, so floating-point types are not recommended.

Algorithm

Similarity transformations first reduce the matrix to upper Hessenberg form. These transformations preserve the characteristic polynomial. The characteristic polynomials of all leading principal submatrices are then built using the Hessenberg recurrence.

Example

#include "math/matrix/characteristic_polynomial.hpp"
#include "math/modint.hpp"

#include <iostream>
#include <vector>

int main() {
    using mint = m1une::math::modint998244353;
    std::vector<mint> values{1, 2, 3, 4};
    m1une::matrix::Matrix<mint> matrix(2, 2, values);

    // det(xI - matrix) = x^2 - 5x - 2.
    std::vector<mint> polynomial =
        m1une::matrix::characteristic_polynomial(matrix);
    for (mint coefficient : polynomial) {
        std::cout << coefficient << ' ';
    }
}

Depends on

Required by

Verified with

Code

#ifndef M1UNE_MATRIX_CHARACTERISTIC_POLYNOMIAL_HPP
#define M1UNE_MATRIX_CHARACTERISTIC_POLYNOMIAL_HPP 1

#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

#include "matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
std::vector<T> characteristic_polynomial(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();

    for (int col = 0; col + 2 < size; col++) {
        int pivot = col + 1;
        while (pivot < size && matrix[pivot][col] == T()) pivot++;
        if (pivot == size) continue;

        if (pivot != col + 1) {
            matrix.swap_rows(pivot, col + 1);
            for (int row = 0; row < size; row++) {
                std::swap(matrix[row][pivot], matrix[row][col + 1]);
            }
        }

        const T inverse_pivot = T(1) / matrix[col + 1][col];
        for (int row = col + 2; row < size; row++) {
            if (matrix[row][col] == T()) continue;
            const T factor = matrix[row][col] * inverse_pivot;
            for (int j = col; j < size; j++) {
                matrix[row][j] -= factor * matrix[col + 1][j];
            }
            for (int i = 0; i < size; i++) {
                matrix[i][col + 1] += factor * matrix[i][row];
            }
        }
    }

    std::vector<std::vector<T>> polynomial(std::size_t(size + 1));
    polynomial[0].assign(1, T(1));
    for (int leading_size = 1; leading_size <= size; leading_size++) {
        const int last = leading_size - 1;
        polynomial[std::size_t(leading_size)].assign(
            std::size_t(leading_size + 1),
            T()
        );
        const std::vector<T>& previous =
            polynomial[std::size_t(leading_size - 1)];
        std::vector<T>& current = polynomial[std::size_t(leading_size)];

        for (int degree = 0; degree < leading_size; degree++) {
            current[std::size_t(degree)] -=
                previous[std::size_t(degree)] * matrix[last][last];
            current[std::size_t(degree + 1)] +=
                previous[std::size_t(degree)];
        }

        T subdiagonal_product = T(1);
        for (int row = last - 1; row >= 0; row--) {
            subdiagonal_product *= matrix[row + 1][row];
            const T factor = subdiagonal_product * matrix[row][last];
            if (factor == T()) continue;
            for (int degree = 0; degree <= row; degree++) {
                current[std::size_t(degree)] -=
                    factor * polynomial[std::size_t(row)][std::size_t(degree)];
            }
        }
    }
    return polynomial[std::size_t(size)];
}

}  // namespace matrix
}  // namespace m1une

#endif  // M1UNE_MATRIX_CHARACTERISTIC_POLYNOMIAL_HPP
#line 1 "math/matrix/characteristic_polynomial.hpp"



#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

#line 1 "math/matrix/matrix.hpp"



#line 6 "math/matrix/matrix.hpp"
#include <cstdint>
#line 9 "math/matrix/matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
class Matrix {
   private:
    int _rows;
    int _cols;
    std::vector<T> _data;

    static std::size_t storage_size(int rows, int cols) {
        assert(rows >= 0);
        assert(cols >= 0);
        return std::size_t(rows) * std::size_t(cols);
    }

   public:
    using value_type = T;

    Matrix() : _rows(0), _cols(0) {}

    Matrix(int rows, int cols, const T& value = T())
        : _rows(rows), _cols(cols), _data(storage_size(rows, cols), value) {}

    Matrix(int rows, int cols, std::vector<T> values)
        : _rows(rows), _cols(cols), _data(std::move(values)) {
        assert(rows >= 0);
        assert(cols >= 0);
        assert(_data.size() == std::size_t(rows) * std::size_t(cols));
    }

    explicit Matrix(const std::vector<std::vector<T>>& values)
        : _rows(int(values.size())), _cols(values.empty() ? 0 : int(values[0].size())),
          _data(storage_size(_rows, _cols)) {
        for (int row = 0; row < _rows; row++) {
            assert(int(values[std::size_t(row)].size()) == _cols);
            for (int col = 0; col < _cols; col++) {
                (*this)[row][col] = values[std::size_t(row)][std::size_t(col)];
            }
        }
    }

    int rows() const {
        return _rows;
    }

    int cols() const {
        return _cols;
    }

    bool empty() const {
        return _rows == 0 || _cols == 0;
    }

    std::vector<T>& data() {
        return _data;
    }

    const std::vector<T>& data() const {
        return _data;
    }

    T* operator[](int row) {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    const T* operator[](int row) const {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    T& operator()(int row, int col) {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    const T& operator()(int row, int col) const {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    static Matrix identity(int size) {
        assert(size >= 0);
        Matrix result(size, size);
        for (int i = 0; i < size; i++) result[i][i] = T(1);
        return result;
    }

    Matrix transposed() const {
        Matrix result(_cols, _rows);
        for (int row = 0; row < _rows; row++) {
            for (int col = 0; col < _cols; col++) {
                result[col][row] = (*this)[row][col];
            }
        }
        return result;
    }

    void swap_rows(int first, int second) {
        assert(0 <= first && first < _rows);
        assert(0 <= second && second < _rows);
        if (first == second) return;
        for (int col = 0; col < _cols; col++) {
            std::swap((*this)[first][col], (*this)[second][col]);
        }
    }

    Matrix& operator+=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] += rhs._data[i];
        return *this;
    }

    Matrix& operator-=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] -= rhs._data[i];
        return *this;
    }

    Matrix& operator*=(const T& scalar) {
        for (T& value : _data) value *= scalar;
        return *this;
    }

    Matrix& operator/=(const T& scalar) {
        for (T& value : _data) value /= scalar;
        return *this;
    }

    Matrix& operator*=(const Matrix& rhs) {
        return *this = *this * rhs;
    }

    Matrix operator+() const {
        return *this;
    }

    Matrix operator-() const {
        Matrix result = *this;
        for (T& value : result._data) value = T() - value;
        return result;
    }

    friend Matrix operator+(Matrix lhs, const Matrix& rhs) {
        return lhs += rhs;
    }

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

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

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

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

    friend Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
        assert(lhs._cols == rhs._rows);
        Matrix result(lhs._rows, rhs._cols);
        for (int row = 0; row < lhs._rows; row++) {
            T* output = result[row];
            for (int middle = 0; middle < lhs._cols; middle++) {
                const T coefficient = lhs[row][middle];
                if (coefficient == T()) continue;
                const T* input = rhs[middle];
                for (int col = 0; col < rhs._cols; col++) {
                    output[col] += coefficient * input[col];
                }
            }
        }
        return result;
    }

    friend std::vector<T> operator*(const Matrix& lhs, const std::vector<T>& rhs) {
        assert(lhs._cols == int(rhs.size()));
        std::vector<T> result(std::size_t(lhs._rows));
        for (int row = 0; row < lhs._rows; row++) {
            T value = T();
            for (int col = 0; col < lhs._cols; col++) {
                value += lhs[row][col] * rhs[std::size_t(col)];
            }
            result[std::size_t(row)] = value;
        }
        return result;
    }

    friend std::vector<T> operator*(const std::vector<T>& lhs, const Matrix& rhs) {
        assert(int(lhs.size()) == rhs._rows);
        std::vector<T> result(std::size_t(rhs._cols));
        for (int row = 0; row < rhs._rows; row++) {
            if (lhs[std::size_t(row)] == T()) continue;
            for (int col = 0; col < rhs._cols; col++) {
                result[std::size_t(col)] += lhs[std::size_t(row)] * rhs[row][col];
            }
        }
        return result;
    }

    bool operator==(const Matrix& rhs) const {
        return _rows == rhs._rows && _cols == rhs._cols && _data == rhs._data;
    }

    bool operator!=(const Matrix& rhs) const {
        return !(*this == rhs);
    }

    Matrix pow(std::uint64_t exponent) const {
        assert(_rows == _cols);
        Matrix result = identity(_rows);
        Matrix base = *this;
        while (exponent > 0) {
            if (exponent & 1) result *= base;
            exponent >>= 1;
            if (exponent > 0) base *= base;
        }
        return result;
    }
};

}  // namespace matrix
}  // namespace m1une


#line 10 "math/matrix/characteristic_polynomial.hpp"

namespace m1une {
namespace matrix {

template <class T>
std::vector<T> characteristic_polynomial(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();

    for (int col = 0; col + 2 < size; col++) {
        int pivot = col + 1;
        while (pivot < size && matrix[pivot][col] == T()) pivot++;
        if (pivot == size) continue;

        if (pivot != col + 1) {
            matrix.swap_rows(pivot, col + 1);
            for (int row = 0; row < size; row++) {
                std::swap(matrix[row][pivot], matrix[row][col + 1]);
            }
        }

        const T inverse_pivot = T(1) / matrix[col + 1][col];
        for (int row = col + 2; row < size; row++) {
            if (matrix[row][col] == T()) continue;
            const T factor = matrix[row][col] * inverse_pivot;
            for (int j = col; j < size; j++) {
                matrix[row][j] -= factor * matrix[col + 1][j];
            }
            for (int i = 0; i < size; i++) {
                matrix[i][col + 1] += factor * matrix[i][row];
            }
        }
    }

    std::vector<std::vector<T>> polynomial(std::size_t(size + 1));
    polynomial[0].assign(1, T(1));
    for (int leading_size = 1; leading_size <= size; leading_size++) {
        const int last = leading_size - 1;
        polynomial[std::size_t(leading_size)].assign(
            std::size_t(leading_size + 1),
            T()
        );
        const std::vector<T>& previous =
            polynomial[std::size_t(leading_size - 1)];
        std::vector<T>& current = polynomial[std::size_t(leading_size)];

        for (int degree = 0; degree < leading_size; degree++) {
            current[std::size_t(degree)] -=
                previous[std::size_t(degree)] * matrix[last][last];
            current[std::size_t(degree + 1)] +=
                previous[std::size_t(degree)];
        }

        T subdiagonal_product = T(1);
        for (int row = last - 1; row >= 0; row--) {
            subdiagonal_product *= matrix[row + 1][row];
            const T factor = subdiagonal_product * matrix[row][last];
            if (factor == T()) continue;
            for (int degree = 0; degree <= row; degree++) {
                current[std::size_t(degree)] -=
                    factor * polynomial[std::size_t(row)][std::size_t(degree)];
            }
        }
    }
    return polynomial[std::size_t(size)];
}

}  // namespace matrix
}  // namespace m1une
Back to top page