la/solve.hpp

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

Depends on

Code

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

// a: row reduced echelon form
template <class T>
vector<vector<T>> ker_basis(const vector<vector<T>>& a, const vector<int>& pivots, int rank) {
  int m = a[0].size();
  vector<vector<T>> res;
  vector<int> ps = {-1};
  rep(i, rank) ps.push_back(pivots[i]);
  ps.push_back(m);
  rep(h, ps.size() - 1) {
    rep2(j, ps[h] + 1, ps[h + 1]) {
      vector<T> x(m);
      x[j] = -T(1);
      rep(i, h) x[pivots[i]] = a[i][j];
      res.push_back(move(x));
    }
  }
  return res;
}

template <class T>
std::optional<pair<vector<T>, vector<vector<T>>>> solve(vector<vector<T>> a, vector<T> b) {
  int n = a.size(), m = a[0].size();
  rep(i, n) a[i].push_back(b[i]);
  auto [pivots, rank] = row_reduce(a);
  rep(i, n) if (pivots[i] == m) return nullopt;
  vector<T> c(n);
  repr(i, rank) c[pivots[i]] = a[i][m];
  return optional(pair(move(c), ker_basis(a, 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);
}
#line 3 "la/solve.hpp"

// a: row reduced echelon form
template <class T>
vector<vector<T>> ker_basis(const vector<vector<T>>& a, const vector<int>& pivots, int rank) {
  int m = a[0].size();
  vector<vector<T>> res;
  vector<int> ps = {-1};
  rep(i, rank) ps.push_back(pivots[i]);
  ps.push_back(m);
  rep(h, ps.size() - 1) {
    rep2(j, ps[h] + 1, ps[h + 1]) {
      vector<T> x(m);
      x[j] = -T(1);
      rep(i, h) x[pivots[i]] = a[i][j];
      res.push_back(move(x));
    }
  }
  return res;
}

template <class T>
std::optional<pair<vector<T>, vector<vector<T>>>> solve(vector<vector<T>> a, vector<T> b) {
  int n = a.size(), m = a[0].size();
  rep(i, n) a[i].push_back(b[i]);
  auto [pivots, rank] = row_reduce(a);
  rep(i, n) if (pivots[i] == m) return nullopt;
  vector<T> c(n);
  repr(i, rank) c[pivots[i]] = a[i][m];
  return optional(pair(move(c), ker_basis(a, pivots, rank)));
}