math/matrix/matrix.hpp

• View this file on GitHub
• Last update: 2024-10-04 14:05:18+02:00
• Include: #include "math/matrix/matrix.hpp"

Depends on

• structure/others/md-vector.hpp

Code

#include <bits/stdc++.h>

#include "../../structure/others/md-vector.hpp"

using namespace std;
template <typename T>
struct Matrix {
    md_vector<T, 2> A;

    Matrix() {}
    Matrix(int n, int m) : A({n, m}) {}
    Matrix(int n) : A({n, n}) {}

    int height() const {
        return A.dsize_[0];
    }

    int width() const {
        return A.dsize_[1];
    }

    inline const auto operator[](int k) const {
        return A[k];
    }
    inline auto operator[](int k) {
        return A[k];
    }

    static Matrix I(int n) {
        Matrix<T> mat(n);
        for (int i = 0; i < n; i++) {
            mat[i][i] = 1;
        }
        return mat;
    }

    Matrix &operator+=(const Matrix &b) {
        int n = height(), m = width();
        assert(n == b.height() && m == b.width());
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                (*this)[i][j] += b[i][j];
            }
        }
        return (*this);
    }

    Matrix &operator-=(Matrix &b) {
        int n = height(), m = width();
        assert(n == b.height() && m == b.width());
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                (*this)[i][j] -= b[i][j];
            }
        }
        return (*this);
    }

    Matrix &operator*=(const Matrix &b) {
        int n = height(), m = b.width(), p = width();

        assert(p == b.height());

        Matrix<T> C(n, m);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                for (int k = 0; k < p; k++) {
                    C[i][j] += (*this)[i][k] * b[k][j];
                }
            }
        }
        swap(this->A, C.A);
        return (*this);
    }

    Matrix &operator^=(long long k) {
        Matrix<T> b = Matrix<T>::I(height());
        while (k) {
            if (k & 1) b *= *this;
            *this *= *this;
            k >>= 1LL;
        }
        swap(this->A, b.A);
        return (*this);
    }

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

    Matrix operator-(const Matrix &b) {
        return (Matrix(*this) -= b);
    }
    Matrix operator*(const Matrix &b) const {
        return (Matrix(*this) *= b);
    }
    Matrix operator^(const long long k) {
        return (Matrix(*this) ^= k);
    }

    T determinant() {
        // 行列式辗转相除法
        assert(this->height() == this->width());
        Matrix<T> B(*this);
        T res = 1;
        int n = this->height();
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                while (B[j][i] != 0) {
                    T t = B[i][i] / B[j][i];
                    for (int k = i; k < n; k++) {
                        B[i][k] -= B[j][k] * t;
                        swap(B[i][k], B[j][k]);
                    }
                    res = -res;
                }
            }
            if (B[i][i] == 0) return 0;
            res *= B[i][i];
        }
        return res;
    }
};

#line 1 "math/matrix/matrix.hpp"
#include <bits/stdc++.h>

#line 3 "structure/others/md-vector.hpp"

#line 5 "structure/others/md-vector.hpp"
using namespace std;

// 多维动态大小数组，可以用于DP等场景。
template <typename T, size_t dimensions>
class md_vector;

namespace internal {
template <size_t dimensions>
size_t md_size(const array<size_t, dimensions>& dsize) {
    size_t base = 1;
    for (int i = 0; i < dimensions; i++) {
        base *= dsize[i];
    }
    return base;
}

template <typename T, size_t dimensions, size_t idx_dimensions>
class md_vector_index {
   public:
    md_vector_index(md_vector<T, dimensions>& vec, size_t base = 0) : vector_(vec), base_(base) {
    }

    auto operator[](size_t v) {
        assert(v < vector_.dsize_[idx_dimensions - 1]);
        return md_vector_index<T, dimensions, idx_dimensions + 1>(vector_, (base_ + v) * vector_.dsize_[idx_dimensions]);
    }

   private:
    md_vector<T, dimensions>& vector_;
    const size_t base_;
};

template <typename T, size_t dimensions>
class md_vector_index<T, dimensions, dimensions> {
   public:
    md_vector_index(md_vector<T, dimensions>& vec, size_t base = 0) : vector_(vec), base_(base) {
    }

    T& operator[](size_t v) {
        return vector_.data_[base_ + v];
    }

    md_vector<T, dimensions>& vector_;
    const size_t base_;
};

template <typename T, size_t dimensions, size_t idx_dimensions>
class const_md_vector_index {
   public:
    const_md_vector_index(const md_vector<T, dimensions>& vec, size_t base = 0) : vector_(vec), base_(base) {
    }

    auto operator[](size_t v) const {
        assert(v < vector_.dsize_[idx_dimensions - 1]);
        return const_md_vector_index<T, dimensions, idx_dimensions + 1>(vector_, (base_ + v) * vector_.dsize_[idx_dimensions]);
    }

   private:
    const md_vector<T, dimensions>& vector_;
    const size_t base_;
};

template <typename T, size_t dimensions>
class const_md_vector_index<T, dimensions, dimensions> {
   public:
    const_md_vector_index(const md_vector<T, dimensions>& vec, size_t base = 0) : vector_(vec), base_(base) {
    }

    const T& operator[](size_t v) const {
        return vector_.data_[base_ + v];
    }

    const md_vector<T, dimensions>& vector_;
    const size_t base_;
};
}  // namespace internal

template <typename T, size_t dimensions>
class md_vector {
   public:
    md_vector() {}
    md_vector(md_vector<T, dimensions>&& other) : data_(other.data_), dsize_(other.dsize_) {
    }
    md_vector(const md_vector<T, dimensions>& other) : data_(other.data_), dsize_(other.dsize_) {
    }

    md_vector(array<size_t, dimensions> dsize, T default_value = T())
        : dsize_(dsize), data_(internal::md_size(dsize), default_value) {
    }

    md_vector& operator=(md_vector<T, dimensions>&& other) {
        data_ = other.data_;
        dsize_ = other.dsize_;
        return *this;
    }
    md_vector& operator=(const md_vector<T, dimensions>& other) {
        data_ = other.data_;
        dsize_ = other.dsize_;
        return *this;
    }

    auto operator[](size_t v) {
        return internal::md_vector_index<T, dimensions, 1>(*this)[v];
    }

    const auto operator[](size_t v) const {
        return internal::const_md_vector_index<T, dimensions, 1>(*this)[v];
    }

    T& operator[](array<size_t, dimensions> idx) {
        size_t base = 0;
        for (int i = 0; i < dimensions; i++) {
            base *= dsize_[i];
            base += idx[i];
        }
        return data_[base];
    }

    vector<T> data_;
    array<size_t, dimensions> dsize_;
};

template <typename T, size_t dimensions>
istream& operator>>(istream& in, md_vector<T, dimensions>& vec) {
    return in >> vec.data_;
}

template <typename T, size_t dimensions>
void make_md_presum(md_vector<T, dimensions>& vec) {
    size_t diff = 1, base = 0;
    for (int currD = dimensions - 1; currD >= 0; currD--) {
        base = diff * vec.dsize_[currD];
        for (size_t i = 0; i + diff < vec.data_.size(); i++) {
            if (i % base + diff < base) {
                vec.data_[i + diff] += vec.data_[i];
            }
        }
        diff = base;
    }
}

template <typename T, int dimensions, int idx_dimensions>
string to_string(internal::md_vector_index<T, dimensions, idx_dimensions>& vec) {
    int sz = vec.vector_.dsize_[idx_dimensions - 1];
    string s = "{";
    for (int i = 0; i < sz; i++) {
        s += to_string(vec[i]);
        if (i != sz - 1) {
            s += ", ";
        }
    }
    s += "}";
    return s;
}

template <typename T, int dimensions>
string to_string(md_vector<T, dimensions>& vec) {
    int sz = vec.dsize_[0];
    string s = "{";
    for (int i = 0; i < sz; i++) {
        auto it = vec[i];
        s += to_string(it);
        if (i != sz - 1) {
            s += ", ";
        }
    }
    s += "}";
    return s;
}
#line 4 "math/matrix/matrix.hpp"

using namespace std;
template <typename T>
struct Matrix {
    md_vector<T, 2> A;

    Matrix() {}
    Matrix(int n, int m) : A({n, m}) {}
    Matrix(int n) : A({n, n}) {}

    int height() const {
        return A.dsize_[0];
    }

    int width() const {
        return A.dsize_[1];
    }

    inline const auto operator[](int k) const {
        return A[k];
    }
    inline auto operator[](int k) {
        return A[k];
    }

    static Matrix I(int n) {
        Matrix<T> mat(n);
        for (int i = 0; i < n; i++) {
            mat[i][i] = 1;
        }
        return mat;
    }

    Matrix &operator+=(const Matrix &b) {
        int n = height(), m = width();
        assert(n == b.height() && m == b.width());
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                (*this)[i][j] += b[i][j];
            }
        }
        return (*this);
    }

    Matrix &operator-=(Matrix &b) {
        int n = height(), m = width();
        assert(n == b.height() && m == b.width());
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                (*this)[i][j] -= b[i][j];
            }
        }
        return (*this);
    }

    Matrix &operator*=(const Matrix &b) {
        int n = height(), m = b.width(), p = width();

        assert(p == b.height());

        Matrix<T> C(n, m);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                for (int k = 0; k < p; k++) {
                    C[i][j] += (*this)[i][k] * b[k][j];
                }
            }
        }
        swap(this->A, C.A);
        return (*this);
    }

    Matrix &operator^=(long long k) {
        Matrix<T> b = Matrix<T>::I(height());
        while (k) {
            if (k & 1) b *= *this;
            *this *= *this;
            k >>= 1LL;
        }
        swap(this->A, b.A);
        return (*this);
    }

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

    Matrix operator-(const Matrix &b) {
        return (Matrix(*this) -= b);
    }
    Matrix operator*(const Matrix &b) const {
        return (Matrix(*this) *= b);
    }
    Matrix operator^(const long long k) {
        return (Matrix(*this) ^= k);
    }

    T determinant() {
        // 行列式辗转相除法
        assert(this->height() == this->width());
        Matrix<T> B(*this);
        T res = 1;
        int n = this->height();
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                while (B[j][i] != 0) {
                    T t = B[i][i] / B[j][i];
                    for (int k = i; k < n; k++) {
                        B[i][k] -= B[j][k] * t;
                        swap(B[i][k], B[j][k]);
                    }
                    res = -res;
                }
            }
            if (B[i][i] == 0) return 0;
            res *= B[i][i];
        }
        return res;
    }
};

Back to top page