matrix/det_arbitrary_mod.hpp

• View this file on GitHub
• Last update: 2024-04-20 12:21:10+09:00
• Include: #include "matrix/det_arbitrary_mod.hpp"

Depends on

• matrix/base.hpp
• modint/base.hpp

Verified with

• test/matrix/Determinant_of_Matrix_Arbitrary_Mod.test.cpp

Code

#pragma once

#include <cassert>

#include "../matrix/base.hpp"
#include "../modint/base.hpp"

namespace ebi {

template <Modint mint> mint det_arbitrary_mod(matrix<mint> a) {
    assert(a.row_size() == a.column_size());
    bool sgn = true;
    int n = a.row_size();
    for (int r = 0; r < n; r++) {
        if (a[r][r] == 0) {
            for (int i = r + 1; i < n; i++) {
                if (a[i][r] != 0) {
                    sgn = !sgn;
                    a.swap(r, i);
                    break;
                }
            }
        }
        if (a[r][r] == 0) return 0;
        for (int i = r + 1; i < n; i++) {
            while (a[i][r] != 0) {
                long long q = a[r][r].val() / a[i][r].val();
                for (int j = r; j < n; j++) {
                    a[r][j] -= a[i][j] * q;
                }
                a.swap(r, i);
                sgn = !sgn;
            }
        }
    }
    mint d = sgn ? 1 : -1;
    for (int i = 0; i < n; i++) d *= a[i][i];
    return d;
}

}  // namespace ebi

#line 2 "matrix/det_arbitrary_mod.hpp"

#include <cassert>

#line 2 "matrix/base.hpp"

#include <algorithm>
#line 5 "matrix/base.hpp"
#include <iostream>
#include <ranges>
#include <vector>

namespace ebi {

template <class T> struct matrix;

template <class T> matrix<T> identify_matrix(int n) {
    matrix<T> a(n, n);
    for (int i = 0; i < n; i++) {
        a[i][i] = 1;
    }
    return a;
}

template <class T> struct matrix {
  private:
    using Self = matrix<T>;

  public:
    matrix(int n_, int m_, T init_val = 0)
        : n(n_), m(m_), data(n * m, init_val) {}

    matrix(const std::vector<std::vector<T>> &a)
        : n((int)a.size()), m((int)a[0].size()) {
        data = std::vector(n * m);
        for (int i = 0; i < n; i++) {
            std::copy(a[i].begin(), a[i].end(), data.begin() + i * m);
        }
    }

    Self operator+(Self &rhs) const noexcept {
        return Self(*this) += rhs;
    }

    Self operator-(Self &rhs) const noexcept {
        return Self(*this) -= rhs;
    }

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

    Self operator/(Self &rhs) const noexcept {
        return Self(*this) /= rhs;
    }

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

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

    Self &operator+=(Self &rhs) noexcept {
        assert(this->size() == rhs.size());
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                data[i][j] += rhs[i][j];
            }
        }
        return *this;
    }

    Self &operator-=(Self &rhs) noexcept {
        assert(this->size() == rhs.size());
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                data[i][j] -= rhs[i][j];
            }
        }
        return *this;
    }

    Self &operator*=(Self &rhs) noexcept {
        int h = n, w = rhs.column_size();
        assert(m == rhs.row_size());
        Self ret(h, w);
        for (int i = 0; i < h; ++i) {
            for (int k = 0; k < m; ++k) {
                for (int j = 0; j < w; ++j) {
                    ret[i][j] += (*this)[i][k] * rhs[k][j];
                }
            }
        }
        return *this = ret;
    }

    Self &operator/=(const Self &rhs) noexcept {
        auto ret = rhs.inv();
        assert(ret);
        return *this *= ret.value();
    }

    Self &operator*=(const T &rhs) noexcept {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                (*this)[i][j] *= rhs;
            }
        }
        return *this;
    }

    const auto operator[](int i) const {
        return std::ranges::subrange(data.begin() + i * m,
                                     data.begin() + (i + 1) * m);
    }

    auto operator[](int i) {
        return std::ranges::subrange(data.begin() + i * m,
                                     data.begin() + (i + 1) * m);
    }

    void swap(int i, int j) {
        std::swap_ranges(data.begin() + i * m, data.begin() + (i + 1) * m,
                         data.begin() + j * m);
    }

    Self transposition() const {
        Self res(m, n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                res[j][i] = data[i][j];
            }
        }
        return res;
    }

    std::optional<Self> inv() const {
        assert(row_size() == column_size());
        Self a = *this;
        Self b = identify_matrix<T>(n);
        for (int r = 0; r < n; r++) {
            for (int i = r; i < n; i++) {
                if (a[i][r] != 0) {
                    a.swap(r, i);
                    b.swap(r, i);
                    break;
                }
            }
            if (a[r][r] == 0) return std::nullopt;
            T x = a[r][r].inv();
            for (int j = 0; j < n; j++) {
                if (r < j) a[r][j] *= x;
                b[r][j] *= x;
            }
            for (int i = 0; i < n; i++) {
                if (i == r) continue;
                for (int j = 0; j < n; j++) {
                    if (r < j) a[i][j] -= a[i][r] * a[r][j];
                    b[i][j] -= a[i][r] * b[r][j];
                }
            }
        }
        return b;
    }

    Self pow(long long k) const {
        assert(row_size() == column_size() && k >= 0);
        Self res = identify_matrix<T>(row_size());
        Self x = *this;
        while (k) {
            if (k & 1) res *= x;
            x *= x;
            k >>= 1;
        }
        return res;
    }

    int row_size() const {
        return n;
    }

    int column_size() const {
        return m;
    }

    std::pair<int, int> size() const {
        return {n, m};
    }

  private:
    int n, m;
    std::vector<T> data;
};

template <class T> T det(matrix<T> a) {
    assert(a.row_size() == a.column_size());
    T d = 1;
    int n = a.row_size();
    for (int r = 0; r < n; r++) {
        if (a[r][r] == 0) {
            for (int i = r + 1; i < n; i++) {
                if (a[i][r] != 0) {
                    a.swap(r, i);
                    d = -d;
                }
            }
        }
        if (a[r][r] == 0) return 0;
        d *= a[r][r];
        T inv = a[r][r].inv();
        for (int i = r + 1; i < n; i++) {
            T x = a[i][r] * inv;
            for (int j = r; j < n; j++) {
                a[i][j] -= x * a[r][j];
            }
        }
    }
    return d;
}

template <class T> std::istream &operator>>(std::istream &os, matrix<T> &a) {
    for (int i = 0; i < a.row_size(); i++) {
        for (int j = 0; j < a.column_size(); j++) {
            os >> a[i][j];
        }
    }
    return os;
}

template <class T>
std::ostream &operator<<(std::ostream &os, const matrix<T> &a) {
    for (int i = 0; i < a.row_size(); i++) {
        for (int j = 0; j < a.column_size(); j++) {
            os << a[i][j];
            if (j < a.column_size() - 1) os << ' ';
        }
        if (i < a.row_size() - 1) os << '\n';
    }
    return os;
}

}  // namespace ebi
#line 2 "modint/base.hpp"

#include <concepts>
#line 5 "modint/base.hpp"
#include <utility>

namespace ebi {

template <class T>
concept Modint = requires(T a, T b) {
    a + b;
    a - b;
    a * b;
    a / b;
    a.inv();
    a.val();
    a.pow(std::declval<long long>());
    T::mod();
};

template <Modint mint> std::istream &operator>>(std::istream &os, mint &a) {
    long long x;
    os >> x;
    a = x;
    return os;
}

template <Modint mint>
std::ostream &operator<<(std::ostream &os, const mint &a) {
    return os << a.val();
}

}  // namespace ebi
#line 7 "matrix/det_arbitrary_mod.hpp"

namespace ebi {

template <Modint mint> mint det_arbitrary_mod(matrix<mint> a) {
    assert(a.row_size() == a.column_size());
    bool sgn = true;
    int n = a.row_size();
    for (int r = 0; r < n; r++) {
        if (a[r][r] == 0) {
            for (int i = r + 1; i < n; i++) {
                if (a[i][r] != 0) {
                    sgn = !sgn;
                    a.swap(r, i);
                    break;
                }
            }
        }
        if (a[r][r] == 0) return 0;
        for (int i = r + 1; i < n; i++) {
            while (a[i][r] != 0) {
                long long q = a[r][r].val() / a[i][r].val();
                for (int j = r; j < n; j++) {
                    a[r][j] -= a[i][j] * q;
                }
                a.swap(r, i);
                sgn = !sgn;
            }
        }
    }
    mint d = sgn ? 1 : -1;
    for (int i = 0; i < n; i++) d *= a[i][i];
    return d;
}

}  // namespace ebi

Back to top page