la/matrix.hpp

View this file on GitHub
Last update: 2024-11-09 00:49:07+00:00
Include: #include "la/matrix.hpp"

Depends on

Code

#pragma once
#include "la/inv.hpp"

template <class T>
class matrix : public vector<vector<T>> {
 public:
  using vector<vector<T>>::vector;
  matrix(int n, int m) : vector<vector<T>>(n, vector<T>(m)){ }

  vector<vector<T>>& as_vec() { return (vector<vector<T>>&)*this; }
  const vector<vector<T>>& as_vec() const { return (vector<vector<T>>&)*this; }
  vector<vector<T>> into_vec() { return move(as_vec()); }
  int n() const { return this->size(); }
  int m() const { return (*this)[0].size(); }

  matrix add(const matrix& a) && {
    assert(n() == a.n() && m() == a.m());
    rep(i, n()) rep(j, m()) (*this)[i][j] += a[i][j];
    return move(*this);
  }
  matrix add(const matrix& a) const& { return matrix(*this).add(a); }
  matrix sub(const matrix& a) && {
    assert(n() == a.n() && m() == a.m());
    rep(i, n()) rep(j, m()) (*this)[i][j] -= a[i][j];
    return move(*this);
  }
  matrix sub(const matrix& a) const& { return matrix(*this).sub(a); }
  matrix mul(const matrix& a) const {
    assert(m() == a.n());
    matrix b(n(), a.m());
    rep(i, n()) rep(k, m()) rep(j, a.m()) b[i][j] += (*this)[i][k] * a[k][j];
    return b;
  }
  matrix inv() && {
    assert(invert(as_vec()));
    return move(*this);
  }
  matrix inv() const& { return matrix(*this).inv(); }
};

#line 2 "prelude.hpp"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using vi = vector<int>;
using vvi = vector<vector<int>>;
using vll = vector<ll>;
using vvll = vector<vector<ll>>;
using vc = vector<char>;
#define rep2(i, m, n) for (auto i = (m); i < (n); i++)
#define rep(i, n) rep2(i, 0, n)
#define repr2(i, m, n) for (auto i = (n); i-- > (m);)
#define repr(i, n) repr2(i, 0, n)
#define all(x) begin(x), end(x)
auto ndvec(int n, auto e) { return vector(n, e); }
auto ndvec(int n, auto ...e) { return vector(n, ndvec(e...)); }
auto comp_key(auto&& f) { return [&](auto&& a, auto&& b) { return f(a) < f(b); }; }
auto& max(const auto& a, const auto& b) { return a < b ? b : a; }
auto& min(const auto& a, const auto& b) { return b < a ? b : a; }
#if __cpp_lib_ranges
namespace R = std::ranges;
namespace V = std::views;
#endif
#line 3 "types.hpp"

template <class It>
using val_t = typename iterator_traits<It>::value_type;
#line 3 "la/row_reduce.hpp"

// Pivots are not 1
// Not echelon form
// Entries below pivots are 0
// Returns pivots and rank
template <class T>
pair<vector<int>, int> row_unreduced(vector<vector<T>>& a) {
  int n = a.size(), m = a[0].size();
  vector<int> pivots(n);
  int i = 0, j = n;
  while (i != j) {
    rep(k, i) if (a[i][pivots[k]] != T(0)) {
      T r = a[i][pivots[k]] / a[k][pivots[k]];
      rep(j, a[i].size()) a[i][j] -= a[k][j] * r;
    }
    int p = 0;
    while (p < m && a[i][p] == T(0)) p++;
    if (p < m) pivots[i++] = p;
    else swap(a[i], a[--j]), pivots[j] = p;
  }
  return make_pair(move(pivots), i);
}

// Returns pivots and rank
template <class T>
pair<vector<int>, int> row_reduce(vector<vector<T>>& a) {
  auto [pivots, rank] = row_unreduced(a);
  rep(i, rank) {
    T r = T(1) / a[i][pivots[i]];
    rep2(j, pivots[i], a[i].size()) a[i][j] *= r;
    rep(i2, i) if (a[i2][pivots[i]] != T(0)) {
      r = a[i2][pivots[i]];
      rep2(j, pivots[i], a[i].size()) a[i2][j] -= a[i][j] * r;
    }
  }
  rep(i, rank) rep(j, i) {
    if (pivots[j] < pivots[j + 1]) break;
    swap(pivots[j], pivots[j + 1]);
    swap(a[j], a[j + 1]);
  }
  return make_pair(pivots, rank);
}
#line 3 "la/inv.hpp"

// Return true on failure
template <class T>
bool invert(vector<vector<T>>& a) {
  int n = a.size();
  assert(a[0].size() == n);
  rep(i, n) a[i].insert(a[i].end(), n, T(0)), a[i][n + i]++;
  int r = row_reduce(a).second;
  if (r != n) return true;
  rep(i, n) a[i].erase(a[i].begin(), a[i].begin() + n);
  return false;
}
#line 3 "la/matrix.hpp"

template <class T>
class matrix : public vector<vector<T>> {
 public:
  using vector<vector<T>>::vector;
  matrix(int n, int m) : vector<vector<T>>(n, vector<T>(m)){ }

  vector<vector<T>>& as_vec() { return (vector<vector<T>>&)*this; }
  const vector<vector<T>>& as_vec() const { return (vector<vector<T>>&)*this; }
  vector<vector<T>> into_vec() { return move(as_vec()); }
  int n() const { return this->size(); }
  int m() const { return (*this)[0].size(); }

  matrix add(const matrix& a) && {
    assert(n() == a.n() && m() == a.m());
    rep(i, n()) rep(j, m()) (*this)[i][j] += a[i][j];
    return move(*this);
  }
  matrix add(const matrix& a) const& { return matrix(*this).add(a); }
  matrix sub(const matrix& a) && {
    assert(n() == a.n() && m() == a.m());
    rep(i, n()) rep(j, m()) (*this)[i][j] -= a[i][j];
    return move(*this);
  }
  matrix sub(const matrix& a) const& { return matrix(*this).sub(a); }
  matrix mul(const matrix& a) const {
    assert(m() == a.n());
    matrix b(n(), a.m());
    rep(i, n()) rep(k, m()) rep(j, a.m()) b[i][j] += (*this)[i][k] * a[k][j];
    return b;
  }
  matrix inv() && {
    assert(invert(as_vec()));
    return move(*this);
  }
  matrix inv() const& { return matrix(*this).inv(); }
};