Language: C++

#include <iostream>
#include <vector>
#include <limits.h>
#include <algorithm>
#include <string>
#include <math.h>
#include <limits.h>
#include <queue>
#include <map>
#include <set>
#include <iomanip>
#include <bitset>
#include <cassert>
#include <random>
#include <functional>
#include <stack>
#include <iomanip>
#include <cassert>
using namespace std;
long long MOD = 1e9 + 7;

//置換群を表すクラス
//V[n]がN(n) = {0,1,...,n-1}に対する置換を表す
class Replacement {
public:
	vector<int> V;
	int N;
	int sgn = 1;   //偶置換->1 奇置換->-1
	Replacement(int n) {
		N = n;
		V.resize(N);
		for (n = 0; n < N; n++)V[n] = n;
	}
	Replacement(vector<int> v) {
		N = v.size();
		V = v;
		vector<int> VV(V.size());
		for (int n = 0; n < V.size(); n++)VV[n] = n;
		int counter = 0;
		for (int i = 0; i < VV.size(); i++) {
			if (V[i] != VV[i]) {
				for (int x = i + 1; x < VV.size(); x++) {
					if (VV[x] != V[i])continue;
					VV[x] = VV[i];
					VV[i] = V[i];
					break;
				}
				counter++;
			}
		}
		sgn = (counter & 1 ? -1 : 1);
	}
	int operator()(int n) {
		return V[n];
	}
	void print() {
		cout << "(";
		for (int n = 0; n < N - 1; n++) cout << V[n] << ",";
		cout << V[N - 1] << ")" << endl;
	}
	//逆元を返す
	Replacement inv() {
		Replacement i(N);
		i.sgn = sgn;
		for (int n = 0; n < N; n++)i.V[V[n]] = n;
		return i;
	}
};

//互換(a,b) \in G(N)を表す(0-index)
class Replacement2 {
public:
	int N, A, B;   //(A,B) (A<B)
	int sgn = -1;   //奇置換より
	Replacement2(int n, int a, int b) {
		N = n;
		A = min(a, b);
		B = max(a, b);
	}
	void print() {
		cout << "(" << A << "," << B << ")" << endl;
	}
};

//AにBを食わせる return A(B(n))
Replacement operator*(Replacement& A, Replacement2& B) {
	Replacement C(A);
	C.V[B.A] = C.V[B.A] ^ C.V[B.B];
	C.V[B.B] = C.V[B.A] ^ C.V[B.B];
	C.V[B.A] = C.V[B.A] ^ C.V[B.B];
	C.sgn *= -1;
	return C;
}

//A(B(n))を返す
Replacement operator*(Replacement2& A, Replacement2& B) {
	Replacement C(A.N);
	C.V[B.A] = C.V[B.A] ^ C.V[B.B];
	C.V[B.B] = C.V[B.A] ^ C.V[B.B];
	C.V[B.A] = C.V[B.A] ^ C.V[B.B];
	C.V[A.A] = C.V[A.A] ^ C.V[A.B];
	C.V[A.B] = C.V[A.A] ^ C.V[A.B];
	C.V[A.A] = C.V[A.A] ^ C.V[A.B];
	return C;
}

//A(B(n))を返す(使う機会があるのかはわからない)
Replacement operator*(Replacement& A, Replacement& B) {
	Replacement C(A.N);
	for (int n = 0; n < C.N; n++) {
		C.V[n] = A.V[B.V[n]];
	}
	C.sgn = A.sgn * B.sgn;
	return C;
}

//ei1333様の行列ライブラリを独自改良させていただきました。ありがとうございます
template< class T >
struct Matrix {
	vector< vector< T > > A;

	Matrix() {}

	Matrix(size_t n, size_t m) : A(n, vector< T >(m, 0)) {}

	Matrix(size_t n) : A(n, vector< T >(n, 0)) {};

	//置換rから行列を表す
	//M = (e_r(0), e_r(1),..., e_r(n-1) )
	//(M(r_1 )*M(r_2 )=M(r_1 (r_2 ))
	Matrix(Replacement r) : A(r.N, vector< T >(r.N, 0)) {
		for (int x = 0; x < A[0].size(); x++) {
			A[r(x)][x] = 1;
		}
	};

	size_t height() const {
		return (A.size());
	}

	size_t width() const {
		return (A[0].size());
	}

	inline const vector< T >& operator[](int k) const {
		return (A.at(k));
	}

	inline vector< T >& operator[](int k) {
		return (A.at(k));
	}

	static Matrix I(size_t n) {
		Matrix mat(n);
		for (int i = 0; i < n; i++) mat[i][i] = 1;
		return (mat);
	}

	Matrix& operator+=(const Matrix& B) {
		size_t n = height(), m = width();
		assert(n == B.height() && m == B.width());
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				(*this)[i][j] += B[i][j];
		return (*this);
	}

	Matrix& operator-=(const Matrix& B) {
		size_t n = height(), m = width();
		assert(n == B.height() && m == B.width());
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				(*this)[i][j] -= B[i][j];
		return (*this);
	}

	Matrix& operator*=(const Matrix& B) {
		size_t n = height(), m = B.width(), p = width();
		assert(p == B.height());
		vector< vector< T > > C(n, vector< T >(m, 0));
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				for (int k = 0; k < p; k++)
					C[i][j] = (C[i][j] + (*this)[i][k] * B[k][j]);
		A.swap(C);
		return (*this);
	}

	Matrix& operator^=(long long k) {
		Matrix B = Matrix::I(height());
		while (k > 0) {
			if (k & 1) B *= *this;
			*this *= *this;
			k >>= 1LL;
		}
		A.swap(B.A);
		return (*this);
	}

	Matrix& operator%=(long long mod) {
		for (int y = 0; y < A.size(); y++) {
			for (int x = 0; x < A[0].size(); x++) {
				A[y][x] %= mod;
			}
		}
		return (*this);
	}

	Matrix operator+(const Matrix& B) const {
		return (Matrix(*this) += B);
	}

	Matrix operator-(const Matrix& B) const {
		return (Matrix(*this) -= B);
	}

	Matrix operator*(const Matrix& B) const {
		return (Matrix(*this) *= B);
	}

	Matrix operator^(const long long k) const {
		return (Matrix(*this) ^= k);
	}

	friend ostream& operator<<(ostream& os, Matrix& p) {
		size_t n = p.height(), m = p.width();
		for (int i = 0; i < n; i++) {
			os << "[";
			for (int j = 0; j < m; j++) {
				os << p[i][j] << (j + 1 == m ? "]\n" : ",");
			}
		}
		return (os);
	}

	Matrix operator%(const long long mod) const {
		return (Matrix(*this) %= mod);
	}

	T determinant() {
		Matrix B(*this);
		assert(width() == height());
		T ret = 1;
		for (int i = 0; i < width(); i++) {
			int idx = -1;
			for (int j = i; j < width(); j++) {
				if (B[j][i] != 0) idx = j;
			}
			if (idx == -1) return (0);
			if (i != idx) {
				ret *= -1;
				swap(B[i], B[idx]);
			}
			ret *= B[i][i];
			T vv = B[i][i];
			for (int j = 0; j < width(); j++) {
				B[i][j] /= vv;
			}
			for (int j = i + 1; j < width(); j++) {
				T a = B[j][i];
				for (int k = 0; k < width(); k++) {
					B[j][k] -= B[i][k] * a;
				}
			}
		}
		return (ret);
	}
};

//L*U*M_R = Aとなる行列L,U,置換Rを求める O(N^3)
//A:分解する行列
//L,U:分解後の行列(2つ併せて格納) 一応Uの対角成分が1
//R:右から掛ける、置換の行列
//不可能(おそらく非正則と同値)ならばfalseを返す
template<class T>
bool LUdecomposition(Matrix<T>& A, Matrix<T>& L, Matrix<T>& U, Replacement& R) {
	int N = A.W;
	U = new Matrix(A);
	L = Matrix<T>::I(A.W);
	for (int q = 0; q < N - 1; q++) {
		int ix = q;
		while (ix < N) {
			if (U[q][ix] != 0)break;
			ix++;
		}
		if (ix == N)return false;
		if (ix != q) {
			for (int y = 0; y < N; y++) {
				T t = U[y][q];
				U[y][q] = U[y][ix];
				U[y][ix] = t;
			}
			R = R * Replacement2(N, q, ix);
		}
		for (int y = q; y < N; y++)L[y][q] = U[y][q];
		for (int x = q; x < N; x++)U[q][x] /= L[q][q];
		for (int y = q + 1; y < N; y++) {
			for (int x = q; x < N; x++) {
				U[y][x] -= L[y][q] * U[q][x];
			}
		}
	}
	L[N - 1][N - 1] = U[N - 1][N - 1];
	U[N - 1][N - 1] = 1;
	return true;
}

//LUx = Bなる方程式を解く O(N^2)
//(LUx = L(Ux) = Ly = B ---> Ux = y)
template<class T>
Matrix<T> solveLUEquation(Matrix<T>& L, Matrix<T>& U, Matrix<T>& B) {
	int N = L.W;
	//Ly = B
	Matrix<T> Y(N, 1);
	Y[0][0] = B[0][0] / L[0][0];
	for (int i = 1; i < N; i++) {
		T temp = B[i][0];
		for (int j = 0; j < i; j++)temp -= L[i][j] * Y[j][0];
		Y[i][0] = temp / L[i][i];
	}
	//Ux = y
	Matrix<T> X(N, 1);
	X[N - 1][0] = Y[N - 1][0];
	for (int j = N - 2; j >= 0; j--) {
		T temp = Y[j][0];
		for (int t = j + 1; t < N; t++)temp -= U[j][t] * X[t][0];
		X[j][0] = temp;
	}
	return X;
}

//LU分解の結果A=L*U*M_Rから逆行列を計算する O(N^3)
template<class T>
Matrix<T> InverseMatrix(Matrix<T>& L, Matrix<T>& U, Replacement& R) {
	//R^-1
	Matrix<T> RR = Matrix<T>(R.inv());
	int N = U.H;
	//U^-1
	Matrix<T> UU(N);
	for (int j = 0; j < N; j++) {
		for (int i = j; i >= 0; i--) {
			UU[i][j] = (i == j ? 1 : 0);
			for (int t = i + 1; t <= j; t++)UU[i][j] -= U[i][t] * UU[t][j];
		}
	}
	//L^-1
	Matrix<T> LL(N);
	for (int j = 0; j < N; j++) {
		for (int i = j; i < N; i++) {
			if (i == j)LL[i][j] = 1;
			for (int t = j; t < i; t++)LL[i][j] -= L[i][t] * LL[t][j];
			LL[i][j] /= L[i][i];
		}
	}
	LL = UU * LL;
	LL = Matrix<T>(R.inv()) * LL;
	return LL;
}

//行列Aから直接逆行列を求める O(N^3)
//非正則の場合Eを返す
template<class T>
Matrix<T> InverseMatrix(Matrix<T>& A) {
	int N = A.H;
	Matrix<T> U(N), L(N);
	Replacement R(N);
	if (!LUdecomposition(A, L, U, R))return Matrix<T>::I(N);
	return InverseMatrix(L, U, R);
}

long long N, K;

Matrix<long long> modpow(Matrix<long long> base, long long p) {
	if (p == 0)return Matrix<long long>::I(base.width());
	else if (p == 1)return base;
	Matrix<long long> ans = modpow(base, p / 2);
	ans *= ans;
	ans %= MOD;
	if (p & 1)ans *= base;
	ans %= MOD;
	return ans;
}

int main() {
	cin >> N >> K;
	Matrix<long long> M(N, N);
	for (long long y = 0; y < N; y++) {
		for (long long x = 0; x < N; x++) {
			long long a;
			cin >> a;
			if(a == 1)M[x][y] = 1;
		}
	}
	Matrix<long long> ans = modpow(M,K);
	long long count = 0;
	for (long long y = 0; y < N; y++) {
		for (long long x = 0; x < N; x++) {
			count = (count + ans[x][y]) % MOD;
		}
	}
	cout << count << endl;
	return 0;
}

Compiler output:
a.cc: In function 'bool LUdecomposition(Matrix<T>&, Matrix<T>&, Matrix<T>&, Replacement&)':
a.cc:293:31: error: no match for 'operator*' (operand types are 'Replacement' and 'Replacement2')
  293 |                         R = R * Replacement2(N, q, ix);
      |                             ~ ^ ~~~~~~~~~~~~~~~~~~~~~~
      |                             |   |
      |                             |   Replacement2
      |                             Replacement
a.cc:86:13: note: candidate: 'Replacement operator*(Replacement&, Replacement2&)' (near match)
   86 | Replacement operator*(Replacement& A, Replacement2& B) {
      |             ^~~~~~~~
a.cc:86:13: note:   conversion of argument 2 would be ill-formed:
a.cc:96:13: note: candidate: 'Replacement operator*(Replacement2&, Replacement2&)'
   96 | Replacement operator*(Replacement2& A, Replacement2& B) {
      |             ^~~~~~~~
a.cc:96:37: note:   no known conversion for argument 1 from 'Replacement' to 'Replacement2&'
   96 | Replacement operator*(Replacement2& A, Replacement2& B) {
      |                       ~~~~~~~~~~~~~~^
a.cc:108:13: note: candidate: 'Replacement operator*(Replacement&, Replacement&)'
  108 | Replacement operator*(Replacement& A, Replacement& B) {
      |             ^~~~~~~~
a.cc:108:52: note:   no known conversion for argument 2 from 'Replacement2' to 'Replacement&'
  108 | Replacement operator*(Replacement& A, Replacement& B) {
      |                                       ~~~~~~~~~~~~~^