test/src/linear_algebra/array_matrix/matrix_det.test.cpp

• View this file on GitHub
• Last update: 2023-05-11 13:21:36+09:00
• Problem: https://judge.yosupo.jp/problem/matrix_det

Depends on

• Array Matrix (library/linear_algebra/array_matrix.hpp)
• Default Operator (library/util/default_operator.hpp)

Code

#define PROBLEM "https://judge.yosupo.jp/problem/matrix_det"

#include <iostream>
#include <atcoder/modint>

#include "library/linear_algebra/array_matrix.hpp"

using mint = atcoder::modint998244353;
using matrix = suisen::SquareArrayMatrix<mint, 500>;

int main() {
    int n;
    std::cin >> n;
    matrix A = matrix::e1();
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            int val;
            std::cin >> val;
            A[i][j] = val;
        }
    }
    std::cout << A.det().val() << '\n';
    return 0;
}

#line 1 "test/src/linear_algebra/array_matrix/matrix_det.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/matrix_det"

#include <iostream>
#include <atcoder/modint>

#line 1 "library/linear_algebra/array_matrix.hpp"



#include <array>
#include <cassert>
#include <optional>

#line 1 "library/util/default_operator.hpp"



namespace suisen {
    namespace default_operator {
        template <typename T>
        auto zero() -> decltype(T { 0 }) { return T { 0 }; }
        template <typename T>
        auto one()  -> decltype(T { 1 }) { return T { 1 }; }
        template <typename T>
        auto add(const T &x, const T &y) -> decltype(x + y) { return x + y; }
        template <typename T>
        auto sub(const T &x, const T &y) -> decltype(x - y) { return x - y; }
        template <typename T>
        auto mul(const T &x, const T &y) -> decltype(x * y) { return x * y; }
        template <typename T>
        auto div(const T &x, const T &y) -> decltype(x / y) { return x / y; }
        template <typename T>
        auto mod(const T &x, const T &y) -> decltype(x % y) { return x % y; }
        template <typename T>
        auto neg(const T &x) -> decltype(-x) { return -x; }
        template <typename T>
        auto inv(const T &x) -> decltype(one<T>() / x)  { return one<T>() / x; }
    } // default_operator
    namespace default_operator_noref {
        template <typename T>
        auto zero() -> decltype(T { 0 }) { return T { 0 }; }
        template <typename T>
        auto one()  -> decltype(T { 1 }) { return T { 1 }; }
        template <typename T>
        auto add(T x, T y) -> decltype(x + y) { return x + y; }
        template <typename T>
        auto sub(T x, T y) -> decltype(x - y) { return x - y; }
        template <typename T>
        auto mul(T x, T y) -> decltype(x * y) { return x * y; }
        template <typename T>
        auto div(T x, T y) -> decltype(x / y) { return x / y; }
        template <typename T>
        auto mod(T x, T y) -> decltype(x % y) { return x % y; }
        template <typename T>
        auto neg(T x) -> decltype(-x) { return -x; }
        template <typename T>
        auto inv(T x) -> decltype(one<T>() / x)  { return one<T>() / x; }
    } // default_operator
} // namespace suisen


#line 9 "library/linear_algebra/array_matrix.hpp"

namespace suisen {
    template <
        typename T, size_t N, size_t M,
        T(*_add)(T, T) = default_operator_noref::add<T>, T(*_neg)(T) = default_operator_noref::neg<T>, T(*_zero)() = default_operator_noref::zero<T>,
        T(*_mul)(T, T) = default_operator_noref::mul<T>, T(*_inv)(T) = default_operator_noref::inv<T>, T(*_one)() = default_operator_noref::one<T>
    >
    struct ArrayMatrix : public std::array<std::array<T, M>, N> {
    private:
        template <typename DummyType = void>
        static constexpr bool is_square_v = N == M;
        template <size_t X, size_t Y>
        using MatrixType = ArrayMatrix<T, X, Y, _add, _neg, _zero, _mul, _inv, _one>;
    public:
        using base_type = std::array<std::array<T, M>, N>;
        using container_type = base_type;
        using row_type = std::array<T, M>;

        using base_type::base_type;
        ArrayMatrix(T diag_val = _zero()) {
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) {
                (*this)[i][j] = (i == j ? diag_val : _zero());
            }
        }
        ArrayMatrix(const container_type& c) : base_type{ c } {}
        ArrayMatrix(const std::initializer_list<row_type>& c) {
            assert(c.size() == N);
            size_t i = 0;
            for (const auto& row : c) {
                for (size_t j = 0; j < M; ++j) (*this)[i][j] = row[j];
                ++i;
            }
        }

        static ArrayMatrix e0() { return ArrayMatrix(_zero()); }
        static MatrixType<M, M> e1() { return MatrixType<M, M>(_one()); }

        int size() const {
            static_assert(is_square_v<>);
            return N;
        }
        std::pair<int, int> shape() const { return { N, M }; }
        int row_size() const { return N; }
        int col_size() const { return M; }

        ArrayMatrix operator+() const { return *this; }
        ArrayMatrix operator-() const {
            ArrayMatrix A;
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) A[i][j] = _neg((*this)[i][j]);
            return A;
        }
        friend ArrayMatrix& operator+=(ArrayMatrix& A, const ArrayMatrix& B) {
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) A[i][j] = _add(A[i][j], B[i][j]);
            return A;
        }
        friend ArrayMatrix& operator-=(ArrayMatrix& A, const ArrayMatrix& B) {
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) A[i][j] = _add(A[i][j], _neg(B[i][j]));
            return A;
        }
        template <size_t K>
        friend MatrixType<N, K>& operator*=(ArrayMatrix& A, const MatrixType<M, K>& B) { return A = A * B; }
        friend ArrayMatrix& operator*=(ArrayMatrix& A, const T& val) {
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) A[i][j] = _mul(A[i][j], val);
            return A;
        }
        friend ArrayMatrix& operator/=(ArrayMatrix& A, const ArrayMatrix& B) { static_assert(is_square_v<>); return A *= *B.inv(); }
        friend ArrayMatrix& operator/=(ArrayMatrix& A, const T& val) { return A *= _inv(val); }

        friend ArrayMatrix operator+(ArrayMatrix A, const ArrayMatrix& B) { A += B; return A; }
        friend ArrayMatrix operator-(ArrayMatrix A, const ArrayMatrix& B) { A -= B; return A; }
        template <size_t K>
        friend MatrixType<N, K> operator*(const ArrayMatrix& A, const MatrixType<M, K>& B) {
            MatrixType<N, K> C;
            for (size_t i = 0; i < N; ++i) {
                C[i].fill(_zero());
                for (size_t j = 0; j < M; ++j) for (size_t k = 0; k < K; ++k) C[i][k] = _add(C[i][k], _mul(A[i][j], B[j][k]));
            }
            return C;
        }
        friend ArrayMatrix operator*(ArrayMatrix A, const T& val) { A *= val; return A; }
        friend ArrayMatrix operator*(const T& val, ArrayMatrix A) {
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) A[i][j] = _mul(val, A[i][j]);
            return A;
        }
        friend std::array<T, N> operator*(const ArrayMatrix& A, const std::array<T, M>& x) {
            std::array<T, N> b;
            b.fill(_zero());
            for (size_t i = 0; i < N; ++i) for (size_t j = 0; j < M; ++j) b[i] = _add(b[i], _mul(A[i][j], x[j]));
            return b;
        }
        friend ArrayMatrix operator/(ArrayMatrix A, const ArrayMatrix& B) { static_assert(is_square_v<>); return A * B.inv(); }
        friend ArrayMatrix operator/(ArrayMatrix A, const T& val) { A /= val; return A; }
        friend ArrayMatrix operator/(const T& val, ArrayMatrix A) { return A.inv() *= val; }

        ArrayMatrix pow(long long b) const {
            static_assert(is_square_v<>);
            assert(b >= 0);
            ArrayMatrix res(e1()), p(*this);
            for (; b; b >>= 1) {
                if (b & 1) res *= p;
                p *= p;
            }
            return res;
        }

        std::optional<ArrayMatrix> safe_inv() const {
            static_assert(is_square_v<>);
            std::array<std::array<T, 2 * N>, N> data;
            for (size_t i = 0; i < N; ++i) {
                for (size_t j = 0; j < N; ++j) {
                    data[i][j] = (*this)[i][j];
                    data[i][N + j] = i == j ? _one() : _zero();
                }
            }
            for (size_t i = 0; i < N; ++i) {
                for (size_t k = i; k < N; ++k) if (data[k][i] != _zero()) {
                    data[i].swap(data[k]);
                    T c = _inv(data[i][i]);
                    for (size_t j = i; j < 2 * N; ++j) data[i][j] = _mul(c, data[i][j]);
                    break;
                }
                if (data[i][i] == _zero()) return std::nullopt;
                for (size_t k = 0; k < N; ++k) if (k != i and data[k][i] != _zero()) {
                    T c = data[k][i];
                    for (size_t j = i; j < 2 * N; ++j) data[k][j] = _add(data[k][j], _neg(_mul(c, data[i][j])));
                }
            }
            ArrayMatrix res;
            for (size_t i = 0; i < N; ++i) std::copy(data[i].begin() + N, data[i].begin() + 2 * N, res[i].begin());
            return res;
        }
        ArrayMatrix inv() const { return *safe_inv(); }
        T det() const {
            static_assert(is_square_v<>);
            ArrayMatrix A = *this;
            bool sgn = false;
            for (size_t j = 0; j < N; ++j) for (size_t i = j + 1; i < N; ++i) if (A[i][j] != _zero()) {
                std::swap(A[j], A[i]);
                T q = _mul(A[i][j], _inv(A[j][j]));
                for (size_t k = j; k < N; ++k) A[i][k] = _add(A[i][k], _neg(_mul(A[j][k], q)));
                sgn = not sgn;
            }
            T res = sgn ? _neg(_one()) : _one();
            for (size_t i = 0; i < N; ++i) res = _mul(res, A[i][i]);
            return res;
        }
        T det_arbitrary_mod() const {
            static_assert(is_square_v<>);
            ArrayMatrix A = *this;
            bool sgn = false;
            for (size_t j = 0; j < N; ++j) for (size_t i = j + 1; i < N; ++i) {
                for (; A[i][j].val(); sgn = not sgn) {
                    std::swap(A[j], A[i]);
                    T q = A[i][j].val() / A[j][j].val();
                    for (size_t k = j; k < N; ++k) A[i][k] -= A[j][k] * q;
                }
            }
            T res = sgn ? -1 : +1;
            for (size_t i = 0; i < N; ++i) res *= A[i][i];
            return res;
        }
    };
    template <
        typename T, size_t N,
        T(*_add)(T, T) = default_operator_noref::add<T>, T(*_neg)(T) = default_operator_noref::neg<T>, T(*_zero)() = default_operator_noref::zero<T>,
        T(*_mul)(T, T) = default_operator_noref::mul<T>, T(*_inv)(T) = default_operator_noref::inv<T>, T(*_one)() = default_operator_noref::one<T>
    >
    using SquareArrayMatrix = ArrayMatrix<T, N, N, _add, _neg, _zero, _mul, _inv, _one>;
} // namespace suisen


#line 7 "test/src/linear_algebra/array_matrix/matrix_det.test.cpp"

using mint = atcoder::modint998244353;
using matrix = suisen::SquareArrayMatrix<mint, 500>;

int main() {
    int n;
    std::cin >> n;
    matrix A = matrix::e1();
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            int val;
            std::cin >> val;
            A[i][j] = val;
        }
    }
    std::cout << A.det().val() << '\n';
    return 0;
}

Back to top page