Submission #4330463

---

Source Code Expand

#include <bits/stdc++.h>
#define int long long int
#define MOD(x) ((x % MOD_N) + MOD_N) % MOD_N
#define FOR(i,a,b) for(int i=(a);i<(b);++i)
#define FORE(i,a,b) for(int i=(a);i<=(b);++i)
#define RFOR(i,a,b) for(int i=(b)-1;i>=(a);--i)
#define RFORE(i,a,b) for(int i=(b);i>=(a);--i)
#define REP(i,n) FOR(i,0,n)
#define RREP(i,n) RFOR(i,0,n)
#define ALL(c) (c).begin(),(c).end()
#define RALL(c) (c).rbegin(),(c).rend()
#define SORT(c) sort(ALL(c))
#define RSORT(c) sort(RALL(c))
#define SZ(c) (int)((c).size())
#define EACH(i,v) for(auto i=v.begin();i!=v.end();++i)
#define REACH(i,v) for(auto i=v.rbegin();i!=v.rend();++i)
#define LB(c,x) distance((c).begin(),lower_bound(ALL(c),x))
#define UB(c,x) distance((c).begin(),upper_bound(ALL(c),x))
#define COUNT(c,x) (lower_bound(ALL(c),x)-upper_bound(ALL(c),x))
#define UNIQUE(c) SORT(c); (c).erase(unique(ALL(c)),(c).end());
#define COPY(c1,c2) copy(ALL(c1),(c2).begin())
#define EXIST(s,e) (bool)((s).find(e)!=(s).end())
#define PB push_back
#define MP make_pair
#define DEL(v) decltype(v)().swap(v)
#define DUMP(x)  cerr<<#x<<" = "<<(x)<<endl;
#define NL cerr<<endl;
using namespace std;
template<typename T,typename U> using P=pair<T,U>;
template<typename T> using V=vector<T>;
template<typename T>bool chmax(T&a,T b){if(a<b){a=b;return true;}return false;}
template<typename T>bool chmin(T&a,T b){if(a>b){a=b;return true;}return false;}
template<typename T>T sum(vector<T>&v){return accumulate(ALL(v),T());}
template<typename T>T sum(vector<T>&v,int a,int b){return accumulate(v.begin()+a,v.begin()+b,T());}
template<typename T>T max(vector<T>&v){return *max_element(ALL(v));}
template<typename T>T min(vector<T>&v){return *max_element(ALL(v));}
template<typename T>T max_index(vector<T>&v){return distance((v).begin(),max_element(ALL(v)));}
template<typename T>T min_index(vector<T>&v){return distance((v).begin(),min_element(ALL(v)));}

struct edge { int to, cost; };

template<typename T>auto&operator<<(ostream&s,const vector<T>&v){s<<"[";bool a=1;for(auto e:v){s<<(a?"":" ")<<e;a=0;}s<<"]";return s;}
template<typename T,typename U>auto&operator<<(ostream&s,const pair<T,U>&p){s<<"("<<p.first<<","<<p.second<<")";return s;}
template<typename T>auto&operator<<(ostream&s,const set<T>&st){s<<"{";bool a=1;for(auto e:st){s<<(a?"":" ")<<e;a=0;}s<<"}";return s;}
template<typename T,typename U>auto&operator<<(ostream&s,const map<T,U>&m){s<<"{";bool a=1;for(auto e:m){s<<(a?"":" ")<<e.first<<":"<<e.second;a=0;}s<<"}";return s;}

const int INF = 1e18;
const int MOD_N = 1e9+7;

template<typename T> class Matrix
{
   private:
      vector<vector<T>> val;
      int m, n;
   public:
      Matrix(int a=0, int b=0) : m(a), n(b) {
         val.resize(m, vector<T>(n));
      }
      Matrix(const vector<vector<T>>& vec) : m(vec.size()), n(vec[0].size()) {
         val.resize(m, vector<T>(n));
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            val[i][j] = vec[i][j];
      }
      ~Matrix() {
         vector<vector<T>>().swap(val);
      }
      Matrix& operator=(const Matrix& mat) {
         try { if (m != mat.m || n != mat.n) throw "Matrix can't be substituted"; }
         catch (const char* e) { cerr << e << endl; }
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            val[i][j] = mat.val[i][j];
         return *this;
      }
      T operator()(int i, int j) {
         try { if (i < 0 || m <= i || j < 0 || n <= j) throw "Bad access"; }
         catch (const char* e) { cerr << e << endl; }
         return val[i][j];
      }
      const Matrix operator+(const Matrix& mat) const {
         try { if (m != mat.m || n != mat.n) throw "Matrix can't be added"; }
         catch (const char* e) { cerr << e << endl; }
         Matrix tmp(m, n);
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            tmp.val[i][j] = val[i][j] + mat.val[i][j];
         return tmp;
      }
      const Matrix operator-(const Matrix& mat) const {
         try { if (m != mat.m || n != mat.n) throw "Matrix can't be subtracted"; }
         catch (const char* e) { cerr << e << endl; }
         Matrix tmp(m, n);
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            tmp.val[i][j] = val[i][j] - mat.val[i][j];
         return tmp;
      }
      const Matrix operator*(const Matrix& mat) const {
         try { if (n != mat.m) throw "Matrix can't be producted"; }
         catch (const char* e) { cerr << e << endl; }
         Matrix tmp(m, mat.n);
         for (int i = 0; i < m; i++) for (int j = 0; j < mat.n; j++) for (int k = 0; k < n; k++)
            // tmp.val[i][j] += val[i][k]*mat.val[k][j];
            tmp.val[i][j] = MOD(tmp.val[i][j] + val[i][k]*mat.val[k][j]);
         return tmp;
      }
      Matrix& operator+=(const Matrix& mat) {
         try { if (m != mat.m || n != mat.n) throw "Matrix can't be added"; }
         catch (const char* e) { cerr << e << endl; }
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            val[i][j] += mat.val[i][j];
         return *this;
      }
      Matrix& operator-=(const Matrix& mat) {
         try { if (m != mat.m || n != mat.n) throw "Matrix can't be subtracted"; }
         catch (const char* e) { cerr << e << endl; }
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            val[i][j] -= mat.val[i][j];
         return *this;
      }
      Matrix operator*=(const Matrix& mat) {
         try { if (n != mat.m) throw "Matrix can't be producted"; }
         catch (const char* e) { cerr << e << endl; }
         Matrix tmp(m, mat.n);
         for (int i = 0; i < m; i++) for (int j = 0; j < mat.n; j++) for (int k = 0; k < n; k++)
            // tmp.val[i][j] += val[i][k]*mat.val[k][j];
            tmp.val[i][j] = MOD(tmp.val[i][j] + val[i][k]*mat.val[k][j]);
         n = mat.n;
         val.resize(m, vector<T>(n));
         *this = tmp;
         return *this;
      }
      Matrix transpose() const {
         Matrix tmp(n, m);
         for (int i = 0; i < m; i++) for (int j = 0; j < n; j++)
            tmp.val[j][i] = val[i][j];
         return tmp;
      }
      Matrix gaussElimination() const {
         Matrix<T> A(val);
         int col = 0;
         for (int i = 0; i < m; i++) {
            int pivot = i;
            while (true) {
               for (int j = i+1; j < m; j++) if (abs(A.val[j][col]) > abs(A.val[pivot][col])) pivot = j;
               if (A.val[pivot][col] == 0) {
                  col++;
                  if (col == n) return A;
               } else break;
            }
            swap(A.val[i], A.val[pivot]);
            for (int j = i+1; j < m; j++) {
               T ratio = A.val[j][col] / A.val[i][col];
               for (int k = col; k < n; k++) A.val[j][k] -= A.val[i][k] * ratio;
            }
            col++;
         }
         return A;
      }
      int rank() const {
         Matrix A = gaussElimination();
         int j = 0;
         for (int i = 0; i < m; i++) {
            while (A.val[i][j] == 0) {
               j++; if (j == n) return i;
            }
         }
         return m;
      }
      ostream& info(ostream& s) const {
         s << endl << "[";
         for (int i = 0; i < m; i++) {
            s << "[";
            for (int j = 0; j < n; j++) {
               s << val[i][j];
               if (j < n-1) s << " ";
            }
            s << "]"; if (i < m-1) s << endl;
         }
         s << "]";
         return s;
      }
};
template<typename T> ostream& operator<<(ostream& s, const Matrix<T>& mat) {
   mat.info(s);
   return s;
}


template<typename T>
Matrix<T> power(Matrix<T> A, int n, int h) {
   V<V<int>> E(n, V<int>(n));
   REP(i, n) {
      E[i][i] = 1;
   }
   Matrix<int> res(E);
   while (h > 0) {
      if ((h & 1) == 1) res *= A;
      A *= A;
      h >>= 1;
   }
   return res;
}

signed main()
{
   int H, R; cin >> H >> R;
   V<V<int>> g(R, V<int>(R));
   REP(i, R) {
      REP(j, R) {
         cin >> g[i][j];
      }
   }

   V<V<V<int>>> dp(R, V<V<int>>(1<<R, V<int>(R)));
   REP(k, R) {
      dp[k][1<<k][k] = 1;
      REP(S, 1<<R) {
         REP(i, R) {
            if (!(S >> i & 1)) continue;
            REP(j, R) {
               if (!(S >> j & 1)) continue;
               // dp[k][S][i] += dp[k][S & ~(1<<i)][j] * g[j][i];
               dp[k][S][i] = MOD(dp[k][S][i] + dp[k][S & ~(1<<i)][j] * g[j][i]);
            }
         }
      }
   }
   V<V<int>> num(R, V<int>(R));
   REP(i, R) {
      REP(j, R) {
         REP(S, 1<<R) {
            // num[i][j] += dp[i][S][j];
            num[i][j] = MOD(num[i][j] + dp[i][S][j]);
         }
      }
   }
   Matrix<int> A(num);
   A = power(A, R, H);
   cout << A(0,0) << endl;

   return 0;
}

Submission Info

Judge Result