la/row_reduce.hpp

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

Depends on

Required by

Code

#pragma once
#include "types.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 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);
}