verify/verify-unit-test/inverse-matrix.test.cpp

• View this file on GitHub
• Last update: 2023-05-29 20:50:32+09:00
• Problem: https://judge.yosupo.jp/problem/aplusb

Depends on

• matrix/gauss-elimination.hpp
• matrix/inverse-matrix.hpp
• modint/montgomery-modint.hpp

Code

#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
//
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>
using namespace std;

#include "../../matrix/gauss-elimination.hpp"
#include "../../modint/montgomery-modint.hpp"
#include "../../matrix/inverse-matrix.hpp"

template <typename mint>
void verify(const vector<vector<mint>>& a, const vector<vector<mint>>& b) {
  int N = a.size();
  auto _a = a;
  auto [rank, det] = GaussElimination(_a);

  if (det == 0 or rank != N) {
    assert(det == 0 and rank != N);
    if (!b.empty()) {
      assert(false and "b is not empty is spite of det == 0.");
    }
    return;
  }
  assert(!b.empty());
  vector<vector<mint>> c(N, vector<mint>(N));
  for (int i = 0; i < N; i++) {
    for (int k = 0; k < N; k++) {
      for (int j = 0; j < N; j++) {
        c[i][j] += a[i][k] * b[k][j];
      }
    }
  }

  for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
      assert(c[i][j] == (i == j));
    }
  }
}

using mint = LazyMontgomeryModInt<998244353>;

void test() {
  int testcase = 1000;

  long long _x = 1;
  auto rng = [&]() { return _x = _x * 100000 % 998244353; };

  for (int t = 0; t < testcase; t++) {
    int N = t % 30 + 1;
    vector<vector<mint>> a(N, vector<mint>(N));
    for (auto& v : a) {
      for (auto& x : v) {
        x = rng();
      }
    }
    auto b = inverse_matrix(a);
    verify(a, b);
  }
}

int main() {
  test();

  int a, b;
  cin >> a >> b;
  cout << a + b << endl;
}

#line 1 "verify/verify-unit-test/inverse-matrix.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
//
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>
using namespace std;

#line 2 "matrix/gauss-elimination.hpp"

#include <utility>
#line 5 "matrix/gauss-elimination.hpp"
using namespace std;

// {rank, det(非正方行列の場合は未定義)} を返す
// 型が double や Rational でも動くはず？(未検証)
//
// pivot 候補 : [0, pivot_end)
template <typename T>
std::pair<int, T> GaussElimination(vector<vector<T>> &a, int pivot_end = -1,
                                   bool diagonalize = false) {
  int H = a.size(), W = a[0].size(), rank = 0;
  if (pivot_end == -1) pivot_end = W;
  T det = 1;
  for (int j = 0; j < pivot_end; j++) {
    int idx = -1;
    for (int i = rank; i < H; i++) {
      if (a[i][j] != T(0)) {
        idx = i;
        break;
      }
    }
    if (idx == -1) {
      det = 0;
      continue;
    }
    if (rank != idx) det = -det, swap(a[rank], a[idx]);
    det *= a[rank][j];
    if (diagonalize && a[rank][j] != T(1)) {
      T coeff = T(1) / a[rank][j];
      for (int k = j; k < W; k++) a[rank][k] *= coeff;
    }
    int is = diagonalize ? 0 : rank + 1;
    for (int i = is; i < H; i++) {
      if (i == rank) continue;
      if (a[i][j] != T(0)) {
        T coeff = a[i][j] / a[rank][j];
        for (int k = j; k < W; k++) a[i][k] -= a[rank][k] * coeff;
      }
    }
    rank++;
  }
  return make_pair(rank, det);
}
#line 2 "modint/montgomery-modint.hpp"

template <uint32_t mod>
struct LazyMontgomeryModInt {
  using mint = LazyMontgomeryModInt;
  using i32 = int32_t;
  using u32 = uint32_t;
  using u64 = uint64_t;

  static constexpr u32 get_r() {
    u32 ret = mod;
    for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
    return ret;
  }

  static constexpr u32 r = get_r();
  static constexpr u32 n2 = -u64(mod) % mod;
  static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30");
  static_assert((mod & 1) == 1, "invalid, mod % 2 == 0");
  static_assert(r * mod == 1, "this code has bugs.");

  u32 a;

  constexpr LazyMontgomeryModInt() : a(0) {}
  constexpr LazyMontgomeryModInt(const int64_t &b)
      : a(reduce(u64(b % mod + mod) * n2)){};

  static constexpr u32 reduce(const u64 &b) {
    return (b + u64(u32(b) * u32(-r)) * mod) >> 32;
  }

  constexpr mint &operator+=(const mint &b) {
    if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
    return *this;
  }

  constexpr mint &operator-=(const mint &b) {
    if (i32(a -= b.a) < 0) a += 2 * mod;
    return *this;
  }

  constexpr mint &operator*=(const mint &b) {
    a = reduce(u64(a) * b.a);
    return *this;
  }

  constexpr mint &operator/=(const mint &b) {
    *this *= b.inverse();
    return *this;
  }

  constexpr mint operator+(const mint &b) const { return mint(*this) += b; }
  constexpr mint operator-(const mint &b) const { return mint(*this) -= b; }
  constexpr mint operator*(const mint &b) const { return mint(*this) *= b; }
  constexpr mint operator/(const mint &b) const { return mint(*this) /= b; }
  constexpr bool operator==(const mint &b) const {
    return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
  }
  constexpr bool operator!=(const mint &b) const {
    return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
  }
  constexpr mint operator-() const { return mint() - mint(*this); }
  constexpr mint operator+() const { return mint(*this); }

  constexpr mint pow(u64 n) const {
    mint ret(1), mul(*this);
    while (n > 0) {
      if (n & 1) ret *= mul;
      mul *= mul;
      n >>= 1;
    }
    return ret;
  }

  constexpr mint inverse() const {
    int x = get(), y = mod, u = 1, v = 0, t = 0, tmp = 0;
    while (y > 0) {
      t = x / y;
      x -= t * y, u -= t * v;
      tmp = x, x = y, y = tmp;
      tmp = u, u = v, v = tmp;
    }
    return mint{u};
  }

  friend ostream &operator<<(ostream &os, const mint &b) {
    return os << b.get();
  }

  friend istream &operator>>(istream &is, mint &b) {
    int64_t t;
    is >> t;
    b = LazyMontgomeryModInt<mod>(t);
    return (is);
  }

  constexpr u32 get() const {
    u32 ret = reduce(a);
    return ret >= mod ? ret - mod : ret;
  }

  static constexpr u32 get_mod() { return mod; }
};
#line 2 "matrix/inverse-matrix.hpp"

#line 4 "matrix/inverse-matrix.hpp"

template <typename mint>
vector<vector<mint>> inverse_matrix(const vector<vector<mint>>& a) {
  int N = a.size();
  assert(N > 0);
  assert(N == (int)a[0].size());

  vector<vector<mint>> m(N, vector<mint>(2 * N));
  for (int i = 0; i < N; i++) {
    copy(begin(a[i]), end(a[i]), begin(m[i]));
    m[i][N + i] = 1;
  }

  auto [rank, det] = GaussElimination(m, N, true);
  if (rank != N) return {};

  vector<vector<mint>> b(N);
  for (int i = 0; i < N; i++) {
    copy(begin(m[i]) + N, end(m[i]), back_inserter(b[i]));
  }
  return b;
}
#line 12 "verify/verify-unit-test/inverse-matrix.test.cpp"

template <typename mint>
void verify(const vector<vector<mint>>& a, const vector<vector<mint>>& b) {
  int N = a.size();
  auto _a = a;
  auto [rank, det] = GaussElimination(_a);

  if (det == 0 or rank != N) {
    assert(det == 0 and rank != N);
    if (!b.empty()) {
      assert(false and "b is not empty is spite of det == 0.");
    }
    return;
  }
  assert(!b.empty());
  vector<vector<mint>> c(N, vector<mint>(N));
  for (int i = 0; i < N; i++) {
    for (int k = 0; k < N; k++) {
      for (int j = 0; j < N; j++) {
        c[i][j] += a[i][k] * b[k][j];
      }
    }
  }

  for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
      assert(c[i][j] == (i == j));
    }
  }
}

using mint = LazyMontgomeryModInt<998244353>;

void test() {
  int testcase = 1000;

  long long _x = 1;
  auto rng = [&]() { return _x = _x * 100000 % 998244353; };

  for (int t = 0; t < testcase; t++) {
    int N = t % 30 + 1;
    vector<vector<mint>> a(N, vector<mint>(N));
    for (auto& v : a) {
      for (auto& x : v) {
        x = rng();
      }
    }
    auto b = inverse_matrix(a);
    verify(a, b);
  }
}

int main() {
  test();

  int a, b;
  cin >> a >> b;
  cout << a + b << endl;
}

Back to top page