test/verify/math/matrix/yosupo-rank-of-matrix.test.cpp

• View this file on GitHub
• Last update: 2024-06-20 15:39:33+09:00
• Problem: https://judge.yosupo.jp/problem/matrix_rank

Depends on

• math/external_gcd.hpp
• Matrix(行列) (math/matrix/matrix.hpp)
• modint (math/modint.hpp)
• template/template.hpp

Code

#define PROBLEM "https://judge.yosupo.jp/problem/matrix_rank"

#include "template/template.hpp"
#include "math/modint.hpp"
#include "math/matrix/matrix.hpp"

using mint = modint998244353;


int main() {
  ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
  int N, M; cin >> N >> M;
  if(N == 0 || M == 0) {
    cout << 0 << endl;
    return 0;
  }

  Matrix<mint>A(N,M);

  for(int i=0; i<N; i++) for(int j=0; j<M; j++) cin >> A[i][j];
  cout << A.rank() << endl;
}

#line 1 "test/verify/math/matrix/yosupo-rank-of-matrix.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/matrix_rank"

#line 1 "template/template.hpp"
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int uint;
template<class T> inline bool chmax(T& a, const T& b) {if (a<b) {a=b; return true;} return false;}
template<class T> inline bool chmin(T& a, const T& b) {if (b<a) {a=b; return true;} return false;}
const int INTINF = 1000001000;
const int INTMAX = 2147483647;
const ll LLMAX = 9223372036854775807;
const ll LLINF = 1000000000000000000;
#line 1 "math/modint.hpp"



#line 1 "math/external_gcd.hpp"



#line 5 "math/external_gcd.hpp"

// g,x,y
template<typename T>
constexpr std::tuple<T, T, T> extendedGCD(T a, T b) {
    T x0 = 1, y0 = 0, x1 = 0, y1 = 1;
    while (b != 0) {
        T q = a / b;
        T r = a % b;
        a = b;
        b = r;
        
        T xTemp = x0 - q * x1;
        x0 = x1;
        x1 = xTemp;
        
        T yTemp = y0 - q * y1;
        y0 = y1;
        y1 = yTemp;
    }
    return {a, x0, y0};
}

#line 5 "math/modint.hpp"

template<int MOD>
struct static_modint {
    int value;

    constexpr static_modint() : value(0) {}

    constexpr static_modint(long long v) {
        value = int(((v % MOD) + MOD) % MOD);
    }

    constexpr static_modint& operator+=(const static_modint& other) {
        if ((value += other.value) >= MOD) value -= MOD;
        return *this;
    }

    constexpr static_modint& operator-=(const static_modint& other) {
        if ((value -= other.value) < 0) value += MOD;
        return *this;
    }

    constexpr static_modint& operator*=(const static_modint& other) {
        value = int((long long)value * other.value % MOD);
        return *this;
    }

    constexpr static_modint operator+(const static_modint& other) const {
        return static_modint(*this) += other;
    }

    constexpr static_modint operator-(const static_modint& other) const {
        return static_modint(*this) -= other;
    }

    constexpr static_modint operator*(const static_modint& other) const {
        return static_modint(*this) *= other;
    }

    constexpr static_modint pow(long long exp) const {
        static_modint base = *this, res = 1;
        while (exp > 0) {
            if (exp & 1) res *= base;
            base *= base;
            exp >>= 1;
        }
        return res;
    }

    constexpr static_modint inv() const {
        //return pow(MOD - 2);
        int g,x,y;
        tie(g,x,y) = extendedGCD(value, MOD);
        assert(g==1);
        if (x < 0) {
            x += MOD;
        }
        //cerr << g << " " << x << " " << y << " " << value << endl;
        //assert((((long)x*value)%MOD + MOD)%MOD == 1);
        return x;
    }

    constexpr static_modint& operator/=(const static_modint& other) {
        return *this *= other.inv();
    }

    constexpr static_modint operator/(const static_modint& other) const {
        return static_modint(*this) /= other;
    }

    constexpr bool operator!=(const static_modint& other) const {
        return val() != other.val();
    }

    constexpr bool operator==(const static_modint& other) const {
        return val() == other.val();
    }

    int val() const {
      return this->value;
    }

    friend std::ostream& operator<<(std::ostream& os, const static_modint& mi) {
        return os << mi.value;
    }

    friend std::istream& operator>>(std::istream& is, static_modint& mi) {
        long long x;
        is >> x;
        mi = static_modint(x);
        return is;
    }
};

template <int mod>
using modint = static_modint<mod>;
using modint998244353  = modint<998244353>;
using modint1000000007 = modint<1000000007>;


#line 1 "math/matrix/matrix.hpp"



template <class T>
struct Matrix{
private: 
  vector<vector<T>>vec;
  int N, M;
public:

  Matrix(int _N, int _M) : N(_N), M(_M), vec(vector<vector<T>>(_N, vector<T>(_M))) {
    assert(_N >= 0 && _M >= 0); // 0*0の行列を返したいときもある(逆行列なかったときとか)
  }

  Matrix<T> operator*(const Matrix<T>& rhs) const  {
    assert(M == rhs.N);
    Matrix ret(N,rhs.M);
    for (int i=0; i<N; i++) for (int k=0; k<M; k++) for(int j=0; j<rhs.M; j++) {
      ret.vec[i][j] += vec[i][k] * rhs.vec[k][j];
    } 

    return ret;
  }

  Matrix<T> operator^(unsigned long long k) const {
    assert(N == M);
    Matrix<T> ret(N, N);
    for(int i=0; i<N; i++) ret[i][i] = T(1);

    Matrix<T> base = *this;

    while (k > 0) {
      if (k & 1) {
        ret *= base;
      }

      base *= base;
      k >>= 1; 
    }

    return ret;
  }

  vector<T>& operator[](int i) {
    assert(i < N);
    return vec[i];
  }

  Matrix<T>& operator*=(const Matrix<T>& b) { return (*this) = (*this) * b; }
  Matrix<T>& operator^=(const unsigned long long k) { return (*this) = (*this) ^ k; }

  // さすがにrankを知るのに副作用があるのはヤバいので
  int rank() const {
    Matrix A = *this;
    return A.sweep(M);
  }

  // サイズを返す。N,Mをconstにしたいけどconstにすると*=や^=が面倒になるため、N,Mを非constのprivateにすることでなんとかする。
  pair<int,int> size() const {
    return make_pair(N, M);
  }

  // 逆行列を返す。なければ0*0行列を返す(これはGifted infantsのマネだが、0*0を返す嬉しさはいまいちわかっていない。変えるかも。)
  Matrix<T> inverse() const {
    assert(N == M);
    Matrix A(N, 2*N);
    for(int i=0; i<N; i++) for(int j=0; j<N; j++) A[i][j] = vec[i][j];
    for(int i=0; i<N; i++) A[i][N+i] = T(1);
    int rank = A.sweep(N);
    if (rank < N) return Matrix(0,0);

    Matrix<T> ret(N, N);
    for(int i=0; i<N; i++) for(int j=0; j<N; j++) ret[i][j] = A[i][N+j];

    return ret;
  }  

private:
// 0<= j < varなj列目について掃き出して、rankを返す

int sweep(int var) {
  assert(var <= M);
  int rank = 0;

  for(int col=0; col<var; col++) {
    int pivot = -1;
    for(int row=rank; row<N; row++) {
      // これがdoubleとかなら、
      // if ( && chmax(mx, asb(A[row][col])) ) みたいな条件を付けて、できるだけ絶対値の大きいpivotを選ぶようにする
      if (vec[row][col] != T(0)) {
        pivot = row;
        break; //double なら違う
      }
    }

    if (pivot == -1) continue;
    swap(vec[pivot], vec[rank]);

    T inv = T(1) / vec[rank][col];
    // pivotの行の先頭が1になるように行を定数倍して揃える
    for(int col2=0; col2<M; col2++) {
      vec[rank][col2] *= inv;
    }

    for(int row=0; row<N; row++) {
      // doubleなら、 && A[row:[col] > EPSのときのみこの操作をする
      if (row != rank) {
        T fac = vec[row][col];
        for(int col2=0; col2<M; col2++) {
          vec[row][col2] -= vec[rank][col2] * fac;
        }
      }
    }
    rank++;
  }

  return rank;
}


};


#line 6 "test/verify/math/matrix/yosupo-rank-of-matrix.test.cpp"

using mint = modint998244353;


int main() {
  ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
  int N, M; cin >> N >> M;
  if(N == 0 || M == 0) {
    cout << 0 << endl;
    return 0;
  }

  Matrix<mint>A(N,M);

  for(int i=0; i<N; i++) for(int j=0; j<M; j++) cin >> A[i][j];
  cout << A.rank() << endl;
}

Back to top page