matrix.cpp

#include "../template/includes.cpp"
#include "../template/typedef.cpp"

template <typename T> class Vec {
protected:
  using iterator = typename std::vector<T>::iterator;
  using const_iterator = typename std::vector<T>::const_iterator;
  using reference = T &;
  using const_reference = const T &;
  std::vector<T> v;
  template <typename Unop> Vec<T> unop_new(Unop op) const {
    Vec<T> res(v.size());
    transform(begin(v), end(v), res.begin(), op);
    return res;
  }
  template <typename Binop> Vec<T> &binop(const Vec<T> &r, Binop op) {
    transform(r.begin(), r.end(), v.begin(), v.begin(), op);
    return *this;
  }
  template <typename Binop> Vec<T> binop_new(const Vec<T> &r, Binop op) const {
    Vec<T> res(v.size());
    transform(r.begin(), r.end(), v.begin(), res.begin(), op);
    return res;
  }

public:
  Vec(int n) : v(n) {}
  Vec(int n, const T &val) : v(n, val) {}
  Vec(const std::vector<T> &w) : v(w) {}
  int size() const noexcept { return v.size(); }
  const_iterator begin() const noexcept { return v.begin(); }
  const_iterator end() const noexcept { return v.end(); }
  iterator begin() noexcept { return v.begin(); }
  iterator end() noexcept { return v.end(); }
  reference operator[](int i) { return v[i]; }
  const_reference operator[](int i) const { return v[i]; }
  Vec<T> operator-() const {
    return unop_new([](T val) { return -val; });
  };
  Vec<T> &operator+=(const Vec<T> &r) {
    return binop(r, [](T x, T y) { return x + y; });
  }
  Vec<T> &operator-=(const Vec<T> &r) {
    return binop(r, [](T x, T y) { return x - y; });
  }
  Vec<T> operator+(const Vec<T> &r) const {
    return binop_new(r, [](T x, T y) { return x + y; });
  }
  Vec<T> operator-(const Vec<T> &r) const {
    return binop_new(r, [](T x, T y) { return x - y; });
  }
  T dot(const Vec<T> &r) const {
    return inner_product(v.begin(), v.end(), r.begin(), T(0));
  }
  T norm() const { return this->dot(v); }
  void push_back(const T &r) { v.push_back(r); }
  void concat(const Vec<T> &r) { v.insert(v.end(), r.begin(), r.end()); }
};

template <typename T> class Matrix : public Vec<Vec<T>> {
public:
  using Vec<Vec<T>>::Vec;
  Matrix(int n, int m, const T &val) : Vec<Vec<T>>::Vec(n, Vec<T>(m, val)) {}
  int Y() const { return this->size(); }
  int X() const { return (*this)[0].size(); }
  Matrix<T> transpose() const {
    const int row = Y(), col = X();
    Matrix res(col, row);
    for (int j = 0; j < col; ++j) {
      for (int i = 0; i < row; ++i) {
        res[j][i] = (*this)[i][j];
      }
    }
    return res;
  }
  Matrix<T> operator*(const Matrix<T> &r) const {
    Matrix<T> tr = r.transpose();
    const int row = Y(), col = tr.Y();
    assert(X() == tr.X());
    Matrix<T> res(row, col);
    for (int i = 0; i < row; ++i) {
      for (int j = 0; j < col; ++j) {
        res[i][j] = (*this)[i].dot(tr[j]);
      }
    }
    return res;
  }
  Vec<T> operator*(const Vec<T> &r) const {
    const int row = Y(), col = r.Y();
    assert(r.size() == col);
    Vec<T> res(row);
    for (int i = 0; i < row; ++i) {
      res[i] = (*this)[i].dot(r);
    }
    return res;
  }
  Matrix<T> &operator*=(const Matrix<T> &r) { return *this = *this * r; }
  Matrix<T> operator^(ll n) const {
    const int m = Y();
    assert(m == X());
    Matrix<T> A = *this, res(m, m, 0);
    for (int i = 0; i < m; ++i) res[i][i] = 1;
    while (n > 0) {
      if (n % 2) res *= A;
      A = A * A;
      n /= 2;
    }
    return res;
  }
  void concat_right(const Vec<T> &r) {
    const int n = Y();
    assert(n == r.size());
    for (int i = 0; i < n; ++i) {
      (*this)[i].push_back(r[i]);
    }
  }
  void concat_right(const Matrix<T> &r) {
    const int n = Y();
    assert(n == r.Y());
    for (int i = 0; i < n; ++i) {
      (*this)[i].concat(r[i]);
    }
  }
  void concat_below(const Vec<T> &r) {
    assert(Y() == 0 || X() == r.size());
    this->push_back(r);
  }
  void concat_below(const Matrix<T> &r) {
    assert(Y() == 0 || X() == r.X());
    for (Vec<T> i : r) (*this).push_back(i);
  }
  int rank() const {
    Matrix<T> A = *this;
    if (Y() == 0) return 0;
    const int n = Y(), m = X();
    int r = 0;
    for (int i = 0; r < n && i < m; ++i) {
      int pivot = r;
      for (int j = r + 1; j < n; ++j) {
        if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;
      }
      std::swap(A[pivot], A[r]);
      if (is_zero(A[r][i])) continue;
      for (int k = m - 1; k >= i; --k) A[r][k] = A[r][k] / A[r][i];
      for (int j = r + 1; j < n; ++j) {
        for (int k = m - 1; k >= i; --k) {
          A[j][k] -= A[r][k] * A[j][i];
        }
      }
      ++r;
    }
    return r;
  }
  T det() const {
    const int n = Y();
    if (n == 0) return 1;
    assert(Y() == X());
    Matrix<T> A = *this;
    T D = 1;
    for (int i = 0; i < n; ++i) {
      int pivot = i;
      for (int j = i + 1; j < n; ++j) {
        if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;
      }
      std::swap(A[pivot], A[i]);
      D = D * A[i][i] * T(i != pivot ? -1 : 1);
      if (is_zero(A[i][i])) break;
      for (int j = i + 1; j < n; ++j) {
        for (int k = n - 1; k >= i; --k) {
          A[j][k] -= A[i][k] * A[j][i] / A[i][i];
        }
      }
    }
    return D;
  }
};