test/yosupo/matrix/system_of_linear_equations.test.cpp

View this file on GitHub
Last update: 2024-07-11 22:21:57+09:00
Problem: https://judge.yosupo.jp/problem/system_of_linear_equations

Code

#define PROBLEM "https://judge.yosupo.jp/problem/system_of_linear_equations"
#include "include/mtl/matrix.hpp"
#include "include/mtl/modular.hpp"
#include <bits/stdc++.h>
using namespace std;
using mint = Modular998244353;

int main() {
    int n,m; cin>>n>>m;
    Matrix<mint,500,500> A;
    for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) {
        cin>>A[i][j];
    }
    Matrix<mint,500,1> b;
    for (int i = 0; i < n; i++) {
        cin>>b[i][0];
    }
    auto [rank,ans_basis,suc] = A.solve(b, n, m);
    if (!suc) {
        cout << -1 << endl;
        return 0;
    }
    cout << rank << endl;
    for (int i = 0; i < m; i++) {
        cout << ans_basis[0][i] << ' ';
    }
    cout << endl;
    for (int i = 0; i < (int)rank; i++) {
        for (int j = 0; j < m; j++) {
            cout << ans_basis[i+1][j] << ' ';
        }
        cout << endl;
    }
}

#line 1 "test/yosupo/matrix/system_of_linear_equations.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/system_of_linear_equations"
#line 2 "include/mtl/matrix.hpp"
#include <cstddef>
#include <array>
#include <vector>
#include <numeric>
#include <algorithm>
#include <cmath>
#include <tuple>
#include <cassert>
#include <iostream>

template <typename TYPE, unsigned ROW, unsigned COLUMN>
class Matrix {
 public:
  using value_type = TYPE;
  static constexpr unsigned R = ROW;
  static constexpr unsigned C = COLUMN;
  using row_type = std::array<value_type, COLUMN>;
  using matrix_type = std::array<row_type, ROW>;

 private:
  matrix_type mat_;

 public:
  Matrix() : mat_({}) {}
  Matrix(std::initializer_list<row_type> list) {
    size_t i = 0;
    for (auto &l : list) {
      mat_[i] = l;
      i++;
    }
  }

  static constexpr Matrix I() {
    Matrix ret;
    auto d = std::min(R, C);
    for (unsigned i = 0; i < d; i++) {
      ret[i][i] = 1;
    }
    return ret;
  }

  Matrix<value_type, C, R> T() const {
    Matrix<value_type, C, R> ret;
    for (unsigned i = 0; i < R; i++) {
      for (unsigned j = 0; j < C; j++) {
        ret[j][i] = mat_[i][j];
      }
    }
    return ret;
  }

  const row_type &operator[](size_t i) const { return mat_[i]; }
  row_type &operator[](size_t i) { return mat_[i]; }

  bool operator==(const Matrix &r) const {
    for (unsigned i = 0; i < R; i++) {
      for (unsigned j = 0; j < C; j++) {
        if (mat_[i][j] != r[i][j])
          return false;
      }
    }
    return true;
  }
  bool operator!=(const Matrix &r) const { return !(*this == r); }
  Matrix &operator+=(const Matrix &r) {
    for (unsigned i = 0; i < R; i++) {
      for (unsigned j = 0; j < C; j++) {
        mat_[i][j] += r[i][j];
      }
    }
    return *this;
  }
  Matrix operator+(const Matrix &r) const {
    return Matrix(*this) += r;
  }
  Matrix operator-() const {
    Matrix ret;
    for (unsigned i = 0; i < R; i++) {
      for (unsigned j = 0; j < C; j++) {
        ret[i][j] = -mat_[i][j];
      }
    }
    return ret;
  }
  Matrix &operator-=(const Matrix &r) {
    for (unsigned i = 0; i < R; i++) {
      for (unsigned j = 0; j < C; j++) {
        mat_[i][j] = -r[i][j];
      }
    }
    return *this;
  }
  Matrix operator-(const Matrix &r) const {
    return Matrix(*this) -= r;
  }
  friend Matrix operator*(const value_type &l, const Matrix &r) {
    Matrix ret;
    for (unsigned i = 0; i < R; i++) {
      for (unsigned j = 0; j < C; j++) {
        ret[i][j] = l * r[i][j];
      }
    }
    return ret;
  }
  template <unsigned _COL>
  Matrix<value_type, R, _COL>
  operator*(const Matrix<value_type, C, _COL> &r) const {
    Matrix<value_type, R, _COL> ret;
    for (unsigned i = 0; i < R; i++) {
      for (unsigned k = 0; k < C; k++) {
        auto mik = mat_[i][k];
        for (unsigned j = 0; j < _COL; j++) {
          ret[i][j] += mik * r[k][j];
        }
      }
    }
    return ret;
  }
  Matrix &operator*=(const Matrix &r) { return *this = *this * r; }
  Matrix pow(unsigned long long k) const {
    assert(R == C);
    Matrix u = *this;
    Matrix t = Matrix::I();
    while (k) {
      if (k & 1) {
        t *= u;
      }
      u *= u;
      k >>= 1;
    }
    return t;
  }
  friend Matrix pow(const Matrix &mat, unsigned long long k) {
    return mat.pow(k);
  }
  /**
   * Return Determinant of matrix
   * using Gaussian elimination
  */
  value_type det(unsigned n = R) const {
    value_type ret = 1;
    auto x = *this;
    for (unsigned i = 0; i < n; i++) {
      unsigned index = i;
      while (index < n and x[index][i] == value_type(0)) {
        ++index;
      }
      if (index == n) {
        return 0;
      }
      if (index != i) {
        std::swap(x[index], x[i]);
        ret = -ret;
      }
      auto ixii = value_type(1) / x[i][i];
      for (unsigned j = i + 1; j < n; j++) {
        value_type factor = x[j][i] * ixii;
        for (unsigned k = i + 1; k < n; k++) {
          x[j][k] -= factor * x[i][k];
        }
      }
    }
    for (unsigned i = 0; i < n; i++) {
      ret *= x[i][i];
    }
    return ret;
  }
  // Return rank of left matrix
  unsigned gaussian_elimination(unsigned row=R, unsigned l_column=C, unsigned r_column=0, 
                                std::vector<unsigned>* principal = nullptr) {
    auto w = l_column+r_column;
    unsigned rank = 0;
    if (principal) principal->resize(row+1, l_column);
    for (unsigned i = 0; i < l_column and rank < row; i++) {
      unsigned index = row;
      value_type zero = 0;
      for (unsigned j = rank; j < row; j++) 
        if (mat_[j][i] != zero) {
        index = j;
        break;
      }
      if (index == row) 
        continue;
      if (index != rank) 
        std::swap(mat_[index], mat_[rank]);
      if (principal) (*principal)[rank] = i;
      auto ixii = value_type(1) / mat_[rank][i];
      mat_[rank][i] = 1;
      for (unsigned j = i + 1; j < w; j++) 
        mat_[rank][j] *= ixii;
      for (unsigned j = 0; j < row; j++) if (j != rank) {
        value_type factor = mat_[j][i];
        mat_[j][i] = zero;
        for (unsigned k = i + 1; k < w; k++) {
          mat_[j][k] -= factor * mat_[rank][k];
        }
      }
      ++rank;
    }
    return rank;
  }
  using system_type = Matrix<value_type, R, C+1>;
  using system_result_matrix = Matrix<value_type, R+1, C>;
  std::tuple<unsigned, system_result_matrix, bool> 
  solve(Matrix<value_type, R, 1>& b, unsigned row=R, unsigned column=C) const {
    // x = [A|b]
    system_type x;
    for (unsigned i = 0; i < row; i++) {
      for (unsigned j = 0; j < column; j++) {
        x[i][j] = mat_[i][j];
      }
      x[i][column] = b[i][0];
    }
    // x = [I|ans] ideally
    std::vector<unsigned> col;
    auto rank = x.gaussian_elimination(row, column, 1, &col);
    for (unsigned i = rank; i < row; i++) if (x[i][column] != value_type(0)) {
      return std::make_tuple(column-rank, system_result_matrix{}, false);
    }
    // { x(ans) } 
    // { basis  }
    system_result_matrix res;
    for (unsigned i = 0; i < rank; i++)
      res[0][col[i]] = x[i][column];
    unsigned t = 0;
    col[rank] = column;
    for (; t < col[0]; t++)
      res[t+1][t] = 1;
    for (unsigned i = 0; i < rank; i++) {
      for (unsigned c = col[i]+1; c < col[i+1]; c++) {
        for (unsigned j = 0; j <= i; j++) {
          res[t+1][col[j]] = -x[j][c];
        }
        res[t+1][c] = 1;
        t++;
      }
    }
    assert(t == column-rank);
    return std::make_tuple(column-rank, res, true);
  }
  std::pair<Matrix, bool> inv(unsigned n=R) const {
    // x = [A|I]
    Matrix<value_type, R, C*2> x;
    for (unsigned i = 0; i < n; i++) {
      for (unsigned j = 0; j < n; j++) {
        x[i][j] = mat_[i][j];
      }
      x[i][n+i] = 1;
    }
    // x = [I|A^-1]
    if (x.gaussian_elimination(n, n, n) != n) {
      return std::make_pair(Matrix{}, false);
    }
    Matrix ret;
    for (unsigned i = 0; i < n; i++) {
      for (unsigned j = 0; j < n; j++) {
        ret[i][j] = x[i][n+j];
      }
    }
    return std::make_pair(ret, true);
  }

};
#line 2 "include/mtl/bit_manip.hpp"
#include <cstdint>
#line 4 "include/mtl/bit_manip.hpp"
#if __cplusplus >= 202002L
#ifndef MTL_CPP20
#define MTL_CPP20
#endif
#include <bit>
#endif

namespace bm {

/// Count 1s for each 8 bits
inline constexpr uint64_t popcnt_e8(uint64_t x) {
  x = (x & 0x5555555555555555) + ((x>>1) & 0x5555555555555555);
  x = (x & 0x3333333333333333) + ((x>>2) & 0x3333333333333333);
  x = (x & 0x0F0F0F0F0F0F0F0F) + ((x>>4) & 0x0F0F0F0F0F0F0F0F);
  return x;
}
/// Count 1s
inline constexpr unsigned popcnt(uint64_t x) {
#ifdef MTL_CPP20
  return std::popcount(x);
#else
  return (popcnt_e8(x) * 0x0101010101010101) >> 56;
#endif
}
/// Alias to mtl::popcnt(x)
constexpr unsigned popcount(uint64_t x) {
  return popcnt(x);
}
/// Count trailing 0s. s.t. *11011000 -> 3
inline constexpr unsigned ctz(uint64_t x) {
#ifdef MTL_CPP20
  return std::countr_zero(x);
#else
  return popcnt((x & (-x)) - 1);
#endif
}
/// Alias to mtl::ctz(x)
constexpr unsigned countr_zero(uint64_t x) {
  return ctz(x);
}
/// Count trailing 1s. s.t. *11011011 -> 2
inline constexpr unsigned cto(uint64_t x) {
#ifdef MTL_CPP20
  return std::countr_one(x);
#else
  return ctz(~x);
#endif
}
/// Alias to mtl::cto(x)
constexpr unsigned countr_one(uint64_t x) {
  return cto(x);
}
inline constexpr unsigned ctz8(uint8_t x) {
  return x == 0 ? 8 : popcnt_e8((x & (-x)) - 1);
}
/// [00..0](8bit) -> 0, [**..*](not only 0) -> 1
inline constexpr uint8_t summary(uint64_t x) {
  constexpr uint64_t hmask = 0x8080808080808080ull;
  constexpr uint64_t lmask = 0x7F7F7F7F7F7F7F7Full;
  auto a = x & hmask;
  auto b = x & lmask;
  b = hmask - b;
  b = ~b;
  auto c = (a | b) & hmask;
  c *= 0x0002040810204081ull;
  return uint8_t(c >> 56);
}
/// Extract target area of bits
inline constexpr uint64_t bextr(uint64_t x, unsigned start, unsigned len) {
  uint64_t mask = len < 64 ? (1ull<<len)-1 : 0xFFFFFFFFFFFFFFFFull;
  return (x >> start) & mask;
}
/// 00101101 -> 00111111 -count_1s-> 6
inline constexpr unsigned log2p1(uint8_t x) {
  if (x & 0x80)
    return 8;
  uint64_t p = uint64_t(x) * 0x0101010101010101ull;
  p -= 0x8040201008040201ull;
  p = ~p & 0x8080808080808080ull;
  p = (p >> 7) * 0x0101010101010101ull;
  p >>= 56;
  return p;
}
/// 00101100 -mask_mssb-> 00100000 -to_index-> 5
inline constexpr unsigned mssb8(uint8_t x) {
  assert(x != 0);
  return log2p1(x) - 1;
}
/// 00101100 -mask_lssb-> 00000100 -to_index-> 2
inline constexpr unsigned lssb8(uint8_t x) {
  assert(x != 0);
  return popcnt_e8((x & -x) - 1);
}
/// Count leading 0s. 00001011... -> 4
inline constexpr unsigned clz(uint64_t x) {
#ifdef MTL_CPP20
  return std::countl_zero(x);
#else
  if (x == 0)
    return 64;
  auto i = mssb8(summary(x));
  auto j = mssb8(bextr(x, 8 * i, 8));
  return 63 - (8 * i + j);
#endif
}
/// Alias to mtl::clz(x)
constexpr unsigned countl_zero(uint64_t x) {
  return clz(x);
}
/// Count leading 1s. 11110100... -> 4
inline constexpr unsigned clo(uint64_t x) {
#ifdef MTL_CPP20
  return std::countl_one(x);
#else
  return clz(~x);
#endif
}
/// Alias to mtl::clo(x)
constexpr unsigned countl_one(uint64_t x) {
  return clo(x);
}

inline constexpr unsigned clz8(uint8_t x) {
  return x == 0 ? 8 : 7 - mssb8(x);
}
inline constexpr uint64_t bit_reverse(uint64_t x) {
  x = ((x & 0x00000000FFFFFFFF) << 32) | ((x & 0xFFFFFFFF00000000) >> 32);
  x = ((x & 0x0000FFFF0000FFFF) << 16) | ((x & 0xFFFF0000FFFF0000) >> 16);
  x = ((x & 0x00FF00FF00FF00FF) << 8) | ((x & 0xFF00FF00FF00FF00) >> 8);
  x = ((x & 0x0F0F0F0F0F0F0F0F) << 4) | ((x & 0xF0F0F0F0F0F0F0F0) >> 4);
  x = ((x & 0x3333333333333333) << 2) | ((x & 0xCCCCCCCCCCCCCCCC) >> 2);
  x = ((x & 0x5555555555555555) << 1) | ((x & 0xAAAAAAAAAAAAAAAA) >> 1);
  return x;
}

/// Check if x is power of 2. 00100000 -> true, 00100001 -> false
constexpr bool has_single_bit(uint64_t x) noexcept {
#ifdef MTL_CPP20
  return std::has_single_bit(x);
#else
  return x != 0 && (x & (x - 1)) == 0;
#endif
}

/// Bit width needs to represent x. 00110110 -> 6
constexpr int bit_width(uint64_t x) noexcept {
#ifdef MTL_CPP20
  return std::bit_width(x);
#else
  return 64 - clz(x);
#endif
}

/// Ceil power of 2. 00110110 -> 01000000
constexpr uint64_t bit_ceil(uint64_t x) {
#ifdef MTL_CPP20
  return std::bit_ceil(x);
#else
  if (x == 0) return 1;
  return 1ull << bit_width(x - 1);
#endif
}

/// Floor power of 2. 00110110 -> 00100000
constexpr uint64_t bit_floor(uint64_t x) {
#ifdef MTL_CPP20
  return std::bit_floor(x);
#else
  if (x == 0) return 0;
  return 1ull << (bit_width(x) - 1);
#endif
}

} // namespace bm
#line 5 "include/mtl/modular.hpp"

template <int MOD>
class Modular {
 private:
  unsigned int val_;

 public:
  static constexpr unsigned int mod() { return MOD; }
  template<class T>
  static constexpr unsigned int safe_mod(T v) {
    auto x = (long long)(v%(long long)mod());
    if (x < 0) x += mod();
    return (unsigned int) x;
  }

  constexpr Modular() : val_(0) {}
  template<class T,
      std::enable_if_t<
          std::is_integral<T>::value && std::is_unsigned<T>::value
      > * = nullptr>
  constexpr Modular(T v) : val_(v%mod()) {}
  template<class T,
      std::enable_if_t<
          std::is_integral<T>::value && !std::is_unsigned<T>::value
      > * = nullptr>
  constexpr Modular(T v) : val_(safe_mod(v)) {}

  constexpr unsigned int val() const { return val_; }
  constexpr Modular& operator+=(Modular x) {
    val_ += x.val();
    if (val_ >= mod()) val_ -= mod();
    return *this;
  }
  constexpr Modular operator-() const { return {mod() - val_}; }
  constexpr Modular& operator-=(Modular x) {
    val_ += mod() - x.val();
    if (val_ >= mod()) val_ -= mod();
    return *this;
  }
  constexpr Modular& operator*=(Modular x) {
    auto v = (long long) val_ * x.val();
    if (v >= mod()) v %= mod();
    val_ = v;
    return *this;
  }
  constexpr Modular pow(long long p) const {
    assert(p >= 0);
    Modular t = 1;
    Modular u = *this;
    while (p) {
      if (p & 1)
        t *= u;
      u *= u;
      p >>= 1;
    }
    return t;
  }
  friend constexpr Modular pow(Modular x, long long p) {
    return x.pow(p);
  }
  constexpr Modular inv() const { return pow(mod()-2); }
  constexpr Modular& operator/=(Modular x) { return *this *= x.inv(); }
  constexpr Modular operator+(Modular x) const { return Modular(*this) += x; }
  constexpr Modular operator-(Modular x) const { return Modular(*this) -= x; }
  constexpr Modular operator*(Modular x) const { return Modular(*this) *= x; }
  constexpr Modular operator/(Modular x) const { return Modular(*this) /= x; }
  constexpr Modular& operator++() { return *this += 1; }
  constexpr Modular operator++(int) { Modular c = *this; ++(*this); return c; }
  constexpr Modular& operator--() { return *this -= 1; }
  constexpr Modular operator--(int) { Modular c = *this; --(*this); return c; }

  constexpr bool operator==(Modular x) const { return val() == x.val(); }
  constexpr bool operator!=(Modular x) const { return val() != x.val(); }

  constexpr bool is_square() const {
    return pow((mod()-1)/2) == 1;
  }
  /**
   * Return x s.t. x * x = a mod p
   * reference: https://zenn.dev/peria/articles/c6afc72b6b003c
  */
  constexpr Modular sqrt() const {
    if (!is_square()) 
      throw std::runtime_error("not square");
    auto mod_eight = mod() % 8;
    if (mod_eight == 3 || mod_eight == 7) {
      return pow((mod()+1)/4);
    } else if (mod_eight == 5) {
      auto x = pow((mod()+3)/8);
      if (x * x != *this)
        x *= Modular(2).pow((mod()-1)/4);
      return x;
    } else {
      Modular d = 2;
      while (d.is_square())
        d += 1;
      auto t = mod()-1;
      int s = bm::ctz(t);
      t >>= s;
      auto a = pow(t);
      auto D = d.pow(t);
      int m = 0;
      Modular dt = 1;
      Modular du = D;
      for (int i = 0; i < s; i++) {
        if ((a*dt).pow(1u<<(s-1-i)) == -1) {
          m |= 1u << i;
          dt *= du;
        }
        du *= du;
      }
      return pow((t+1)/2) * D.pow(m/2);
    }
  }

  friend std::ostream& operator<<(std::ostream& os, const Modular& x) {
    return os << x.val();
  }
  friend std::istream& operator>>(std::istream& is, Modular& x) {
    return is >> x.val_;
  }

};

using Modular998244353 = Modular<998244353>;
using Modular1000000007 = Modular<(int)1e9+7>;

template<int Id=0>
class DynamicModular {
 private:
  static unsigned int mod_;
  unsigned int val_;

 public:
  static unsigned int mod() { return mod_; }
  static void set_mod(unsigned int m) { mod_ = m; }
  template<class T>
  static constexpr unsigned int safe_mod(T v) {
    auto x = (long long)(v%(long long)mod());
    if (x < 0) x += mod();
    return (unsigned int) x;
  }

  constexpr DynamicModular() : val_(0) {}
  template<class T,
      std::enable_if_t<
          std::is_integral<T>::value && std::is_unsigned<T>::value
      > * = nullptr>
  constexpr DynamicModular(T v) : val_(v%mod()) {}
  template<class T,
      std::enable_if_t<
          std::is_integral<T>::value && !std::is_unsigned<T>::value
      > * = nullptr>
  constexpr DynamicModular(T v) : val_(safe_mod(v)) {}

  constexpr unsigned int val() const { return val_; }
  constexpr DynamicModular& operator+=(DynamicModular x) {
    val_ += x.val();
    if (val_ >= mod()) val_ -= mod();
    return *this;
  }
  constexpr DynamicModular operator-() const { return {mod() - val_}; }
  constexpr DynamicModular& operator-=(DynamicModular x) {
    val_ += mod() - x.val();
    if (val_ >= mod()) val_ -= mod();
    return *this;
  }
  constexpr DynamicModular& operator*=(DynamicModular x) {
    auto v = (long long) val_ * x.val();
    if (v >= mod()) v %= mod();
    val_ = v;
    return *this;
  }
  constexpr DynamicModular pow(long long p) const {
    assert(p >= 0);
    DynamicModular t = 1;
    DynamicModular u = *this;
    while (p) {
      if (p & 1)
        t *= u;
      u *= u;
      p >>= 1;
    }
    return t;
  }
  friend constexpr DynamicModular pow(DynamicModular x, long long p) {
    return x.pow(p);
  }
  // TODO: implement when mod is not prime
  constexpr DynamicModular inv() const { return pow(mod()-2); }
  constexpr DynamicModular& operator/=(DynamicModular x) { return *this *= x.inv(); }
  constexpr DynamicModular operator+(DynamicModular x) const { return DynamicModular(*this) += x; }
  constexpr DynamicModular operator-(DynamicModular x) const { return DynamicModular(*this) -= x; }
  constexpr DynamicModular operator*(DynamicModular x) const { return DynamicModular(*this) *= x; }
  constexpr DynamicModular operator/(DynamicModular x) const { return DynamicModular(*this) /= x; }
  constexpr DynamicModular& operator++() { return *this += 1; }
  constexpr DynamicModular operator++(int) { DynamicModular c = *this; ++(*this); return c; }
  constexpr DynamicModular& operator--() { return *this -= 1; }
  constexpr DynamicModular operator--(int) { DynamicModular c = *this; --(*this); return c; }

  constexpr bool operator==(DynamicModular x) const { return val() == x.val(); }
  constexpr bool operator!=(DynamicModular x) const { return val() != x.val(); }

  constexpr bool is_square() const {
    return val() == 0 or pow((mod()-1)/2) == 1;
  }
  /**
   * Return x s.t. x * x = a mod p
   * reference: https://zenn.dev/peria/articles/c6afc72b6b003c
  */
  constexpr DynamicModular sqrt() const {
    // assert mod is prime
    if (!is_square()) 
      throw std::runtime_error("not square");
    if (val() < 2)
      return val();
    auto mod_eight = mod() % 8;
    if (mod_eight == 3 || mod_eight == 7) {
      return pow((mod()+1)/4);
    } else if (mod_eight == 5) {
      auto x = pow((mod()+3)/8);
      if (x * x != *this)
        x *= DynamicModular(2).pow((mod()-1)/4);
      return x;
    } else {
      DynamicModular d = 2;
      while (d.is_square())
        ++d;
      auto t = mod()-1;
      int s = bm::ctz(t);
      t >>= s;
      auto a = pow(t);
      auto D = d.pow(t);
      int m = 0;
      DynamicModular dt = 1;
      DynamicModular du = D;
      for (int i = 0; i < s; i++) {
        if ((a*dt).pow(1u<<(s-1-i)) == -1) {
          m |= 1u << i;
          dt *= du;
        }
        du *= du;
      }
      return pow((t+1)/2) * D.pow(m/2);
    }
  }

  friend std::ostream& operator<<(std::ostream& os, const DynamicModular& x) {
    return os << x.val();
  }
  friend std::istream& operator>>(std::istream& is, DynamicModular& x) {
    return is >> x.val_;
  }

};
template<int Id>
unsigned int DynamicModular<Id>::mod_;

#line 264 "include/mtl/modular.hpp"

template<class ModInt>
struct ModularUtil {
  static constexpr int mod = ModInt::mod();
  static struct inv_table {
    std::vector<ModInt> tb{0,1};
    inv_table() : tb({0,1}) {}
  } inv_;
  void set_inv(int n) {
    int m = inv_.tb.size();
    if (m > n) return;
    inv_.tb.resize(n+1);
    for (int i = m; i < n+1; i++)
      inv_.tb[i] = -inv_.tb[mod % i] * (mod / i);
  }
  ModInt& inv(int i) {
    set_inv(i);
    return inv_.tb[i];
  }
};
template<class ModInt>
typename ModularUtil<ModInt>::inv_table ModularUtil<ModInt>::inv_;

#line 288 "include/mtl/modular.hpp"

namespace math {

constexpr int mod_pow_constexpr(int x, int p, int m) {
  long long t = 1;
  long long u = x;
  while (p) {
    if (p & 1) {
      t *= u;
      t %= m;
    }
    u *= u;
    u %= m;
    p >>= 1;
  }
  return (int) t;
}

constexpr int primitive_root_constexpr(int m) {
  if (m == 2) return 1;
  if (m == 167772161) return 3;
  if (m == 469762049) return 3;
  if (m == 754974721) return 11;
  if (m == 880803841) return 26;
  if (m == 998244353) return 3;

  std::array<int, 20> divs{};
  int cnt = 0;
  int x = m-1;
  if (x % 2 == 0) {
    divs[cnt++] = 2;
    x >>= bm::ctz(x);
  }
  for (int d = 3; d*d <= x; d += 2) {
    if (x % d == 0) {
      divs[cnt++] = d;
      while (x % d == 0)
        x /= d;
    }
  }
  if (x > 1) divs[cnt++] = x;
  for (int g = 2; g < m; g++) {
    bool ok = true;
    for (int i = 0; i < cnt; i++) {
      if (mod_pow_constexpr(g, (m-1) / divs[i], m) == 1) {
        ok = false;
        break;
      }
    }
    if (ok) return g;
  }
  return -1;
}

template<int m>
constexpr int primitive_root = primitive_root_constexpr(m);

}
#line 4 "test/yosupo/matrix/system_of_linear_equations.test.cpp"
#include <bits/stdc++.h>
using namespace std;
using mint = Modular998244353;

int main() {
    int n,m; cin>>n>>m;
    Matrix<mint,500,500> A;
    for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) {
        cin>>A[i][j];
    }
    Matrix<mint,500,1> b;
    for (int i = 0; i < n; i++) {
        cin>>b[i][0];
    }
    auto [rank,ans_basis,suc] = A.solve(b, n, m);
    if (!suc) {
        cout << -1 << endl;
        return 0;
    }
    cout << rank << endl;
    for (int i = 0; i < m; i++) {
        cout << ans_basis[0][i] << ' ';
    }
    cout << endl;
    for (int i = 0; i < (int)rank; i++) {
        for (int j = 0; j < m; j++) {
            cout << ans_basis[i+1][j] << ' ';
        }
        cout << endl;
    }
}

Test cases

Env	Name	Status	Elapsed	Memory
g++	example_00	AC	9 ms	7 MB
g++	fullrank_00	AC	45 ms	7 MB
g++	fullrank_01	AC	45 ms	7 MB
g++	fullrank_02	AC	24 ms	7 MB
g++	fullrank_no_solution_00	AC	45 ms	7 MB
g++	fullrank_no_solution_01	AC	46 ms	7 MB
g++	fullrank_no_solution_02	AC	24 ms	7 MB
g++	hack_00	AC	30 ms	7 MB
g++	hack_01	AC	24 ms	7 MB
g++	hack_02	AC	27 ms	7 MB
g++	hack_03	AC	11 ms	7 MB
g++	hack_04	AC	17 ms	7 MB
g++	hack_05	AC	26 ms	7 MB
g++	lowrank_00	AC	117 ms	7 MB
g++	lowrank_01	AC	106 ms	7 MB
g++	lowrank_02	AC	66 ms	7 MB
g++	lowrank_03	AC	95 ms	7 MB
g++	lowrank_04	AC	92 ms	7 MB
g++	lowrank_05	AC	105 ms	7 MB
g++	max_00	AC	176 ms	7 MB
g++	max_01	AC	176 ms	7 MB
g++	max_02	AC	176 ms	7 MB
g++	random_00	AC	37 ms	7 MB
g++	random_01	AC	31 ms	7 MB
g++	random_02	AC	29 ms	7 MB
g++	random_03	AC	31 ms	7 MB
g++	random_04	AC	15 ms	7 MB
clang++	example_00	AC	8 ms	7 MB
clang++	fullrank_00	AC	43 ms	7 MB
clang++	fullrank_01	AC	44 ms	7 MB
clang++	fullrank_02	AC	22 ms	7 MB
clang++	fullrank_no_solution_00	AC	42 ms	6 MB
clang++	fullrank_no_solution_01	AC	43 ms	6 MB
clang++	fullrank_no_solution_02	AC	22 ms	6 MB
clang++	hack_00	AC	23 ms	7 MB
clang++	hack_01	AC	20 ms	7 MB
clang++	hack_02	AC	22 ms	7 MB
clang++	hack_03	AC	9 ms	7 MB
clang++	hack_04	AC	13 ms	7 MB
clang++	hack_05	AC	23 ms	7 MB
clang++	lowrank_00	AC	110 ms	6 MB
clang++	lowrank_01	AC	98 ms	7 MB
clang++	lowrank_02	AC	64 ms	6 MB
clang++	lowrank_03	AC	88 ms	7 MB
clang++	lowrank_04	AC	86 ms	6 MB
clang++	lowrank_05	AC	96 ms	7 MB
clang++	max_00	AC	166 ms	7 MB
clang++	max_01	AC	166 ms	7 MB
clang++	max_02	AC	166 ms	7 MB
clang++	random_00	AC	34 ms	7 MB
clang++	random_01	AC	29 ms	7 MB
clang++	random_02	AC	26 ms	7 MB
clang++	random_03	AC	29 ms	7 MB
clang++	random_04	AC	13 ms	7 MB