Problem:
 Score : 100 points    Problem Statement There is a simple directed graph G with N vertices, numbered 1, 2, \ldots, N .  For each i and j ( 1 \leq i, j \leq N ), you are given an integer a_{i, j} that represents whether there is a directed edge from Vertex i to j .
If a_{i, j} = 1 , there is a directed edge from Vertex i to j ; if a_{i, j} = 0 , there is not.  Find the number of different directed paths of length K in G , modulo 10^9 + 7 .
We will also count a path that traverses the same edge multiple times.      Constraints  All values in input are integers.  1 \leq N \leq 50  1 \leq K \leq 10^{18}  a_{i, j} is 0 or 1 .  a_{i, i} = 0         Input Input is given from Standard Input in the following format:  N  K  a_{1, 1}  \ldots  a_{1, N}  :  a_{N, 1}  \ldots  a_{N, N}       Output Print the number of different directed paths of length K in G , modulo 10^9 + 7 .        Sample Input 1 4 2
0 1 0 0
0 0 1 1
0 0 0 1
1 0 0 0      Sample Output 1 6  G is drawn in the figure below:   There are six directed paths of length 2 :   1 → 2 → 3  1 → 2 → 4  2 → 3 → 4  2 → 4 → 1  3 → 4 → 1  4 → 1 → 2        Sample Input 2 3 3
0 1 0
1 0 1
0 0 0      Sample Output 2 3  G is drawn in the figure below:   There are three directed paths of length 3 :   1 → 2 → 1 → 2  2 → 1 → 2 → 1  2 → 1 → 2 → 3        Sample Input 3 6 2
0 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 0 0
0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 0 0 0      Sample Output 3 1  G is drawn in the figure below:   There is one directed path of length 2 :   4 → 5 → 6        Sample Input 4 1 1
0      Sample Output 4 0       Sample Input 5 10 1000000000000000000
0 0 1 1 0 0 0 1 1 0
0 0 0 0 0 1 1 1 0 0
0 1 0 0 0 1 0 1 0 1
1 1 1 0 1 1 0 1 1 0
0 1 1 1 0 1 0 1 1 1
0 0 0 1 0 0 1 0 1 0
0 0 0 1 1 0 0 1 0 1
1 0 0 0 1 0 1 0 0 0
0 0 0 0 0 1 0 0 0 0
1 0 1 1 1 0 1 1 1 0      Sample Output 5 957538352  Be sure to print the count modulo 10^9 + 7 .

Solution (Python):
N,K = map(int,input().split())
G = [list(map(int,input().split())) for _ in range(N)]

MOD = 10**9+7

def dot(a,b):
  res = [[0] * N for _ in range(N)]
  for i in range(N):
    for j in range(N):
      res[i][j] = sum(a[i][k] * b[k][j] for k in range(N)) % MOD
  return res

ans = [[0] * N for _ in range(N)]
for i in range(N):
  ans[i][i] = 1

while(K>0):
  K,x = divmod(K,2)
  if x == 1:
    ans = dot(ans,G)
  G = dot(G,G)
print(sum(sum(a) % MOD for a in ans) % MOD)

Solution (Python):
import sys

sys.setrecursionlimit(10 ** 6)
int1 = lambda x: int(x) - 1
p2D = lambda x: print(*x, sep="\n")
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(rows_number): return [LI() for _ in range(rows_number)]

def main():
    def dot(aa, bb):
        return [[sum(a * b for a, b in zip(ar, bc)) % md for bc in zip(*bb)] for ar in aa]

    md = 10 ** 9 + 7
    n, k = MI()
    aa = LLI(n)
    ans=[[(i==j)*1 for j in range(n)] for i in range(n)]
    while k:
        if k & 1:
            ans = dot(ans, aa)
        aa = dot(aa, aa)
        k >>= 1
    
    print(sum(sum(ar) for ar in ans) % md)

main()

Solution (Python):


import sys

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 18
MOD = 10 ** 9 + 7

def mat_dot(A, B):
    """ 行列の積 """

    
    if not isinstance(A[0], list) and not isinstance(A[0], tuple):
        A = [A]
    if not isinstance(B[0], list) and not isinstance(A[0], tuple):
        B = [[b] for b in B]
    n1 = len(A)
    n2 = len(A[0])
    _ = len(B)
    m2 = len(B[0])
    res=list2d(n1, m2, 0)
    for i in range(n1):
        for j in range(m2):
            for k in range(n2):
                res[i][j] += A[i][k] * B[k][j]
                res[i][j] %= MOD
    return res

def mat_pow(mat, n, k):
    """ n*nの行列matを行列累乗でk乗する """

    res = list2d(n, n, 0)
    for i in range(n):
        res[i][i] = 1
    
    while k > 0:
        if k & 1:
            res = mat_dot(res, mat)
        mat = mat_dot(mat, mat)
        k >>= 1
    return res

N, K = MAP()
G = [None] * N
for i in range(N):
    G[i] = LIST()


mat = mat_pow(G, N, K)


print(sum([sum(row)%MOD for row in mat])%MOD)

Solution (Python):
def power(a, k):
    N = len(a)
    r = [[1 if i==j else 0 for i in range(N)] for j in range(N)]
    while k:
        if k & 1:
            r = mul(r, a)
        k >>= 1
        a = mul(a, a)
    return r
def mul(a, b):
    return [[sum(ay * by for ay, by in zip(ax, bx)) % (10**9 + 7) for bx in zip(*b)] for ax in a]
def main():
    N, K = map(int, input().split())
    l = []
    for _ in range(N):
        l.append(list(map(int, input().split())))
    lt = [i for i in zip(*l)]
    return sum(mul([[1]*N], power(lt, K))[0])
    

print(main() % (10**9 + 7))

Solution (Python):
def power(a, k):
    N = len(a)
    r = [[1 if i==j else 0 for i in range(N)] for j in range(N)]
    while k:
        if k & 1:
            r = mul(r, a)
        k >>= 1
        a = mul(a, a)
    return r
def mul(a, b):
    return [[sum(ay * by for ay, by in zip(ax, bx)) % (10**9 + 7) for bx in zip(*b)] for ax in a]
def main():
    N, K = map(int, input().split())
    l = []
    for _ in range(N):
        l.append(list(map(int, input().split())))
    lt = [i for i in zip(*l)]
    return sum(mul([[1]*N], power(lt, K))[0])
    

print(main() % (10**9 + 7))

Solution (Python):
MOD = 10**9 + 7

def mul(a, b):
    return [
        [sum(aij * bij for aij, bij in zip(ai, bi)) % MOD for bi in zip(*b)]
        for ai in a]

def pow(b, k):
    r = [[1 if i == j else 0 for i in range(len(b))] for j in range(len(b))]
    while k:
        if k & 1:
            r = mul(r, b)
        k >>= 1
        b = mul(b, b)
    return r


def main():
    N, K = (int(i) for i in input().split())
    A = []
    for i in range(N):
        A.append([int(j) for j in input().split()])
    A = [i for i in zip(*A)]
    v = [[1] * N]
    return sum(mul(v, pow(A, K))[0]) % MOD

print(main())

Solution (Python):
import numpy as np
mod = 10**9+7
def MatPower(A, k):
    B = np.identity(len(A), dtype = object)
    while k:
        if k&1:
            B = np.dot(A, B)%mod
        A = np.dot(A, A)%mod
        k//=2
    return B
N, K = map(int, input().split())
a = np.array([list(map(int, input().split())) for _ in range(N)], dtype = object)
v = np.array([1]*N, dtype = object).T
An = MatPower(a, K)
v = np.dot(An, v)%mod
print(v.sum()%mod)

Solution (Python):
import numpy as np
mod = 10**9+7
def MatPower(A, k):
    B = np.identity(len(A), dtype = object)
    while k:
        if k&1:
            B = np.dot(A, B)%mod
        A = np.dot(A, A)%mod
        k//=2
    return B
N, K = map(int, input().split())
a = np.array([list(map(int, input().split())) for _ in range(N)], dtype = object)
v = np.array([1]*N, dtype = object).T
An = MatPower(a, K)
v = np.dot(An, v)%mod
ans = 0
for i in range(N):
    ans+=v[i]
    ans%=mod
print(ans)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

MOD = 10**9 + 7

N,K = map(int,readline().split())
A = np.array(read().split(),np.int64).reshape(N,N)

def mult(A,B):
    C = np.zeros((N,N),np.int64)
    for n in range(N):
        C[n] = (A[n,:][:,None] * B % MOD).sum(axis=0) % MOD
    return C

def mat_power(A,n):
    if n == 0:
        return np.eye(N,dtype=np.int64)
    X = mat_power(A,n//2)
    X = mult(X,X)
    return mult(A,X) if n&1 else X

answer = mat_power(A,K).sum() % MOD
print(answer)

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=10**9+7
input=lambda :sys.stdin.readline().rstrip()
def power(b,p,f):
    """
    b: object (base)
    p: positive int (multiplier)
    f: multiple function
    return: f^p(b)
    """
    if not isinstance(p,int): raise ValueError("multiplier must be int")
    elif p<=0: raise ValueError("multiplier must be positive.")
    logp=p.bit_length()
    S=[0]*logp
    S[0]=b
    res='$'
    for i in range(logp):
        if i!=logp-1: S[i+1]=f(S[i],S[i])
        if p&(1<<i):
            if res=='$': res=S[i]
            else: res=f(res,S[i])
    return res

def mat_mul(A,B):
    assert len(A[0])==len(B)
    m=len(A)
    n=len(B)
    l=len(B[0])
    C=[[0]*l for _ in range(m)]
    for i in range(m):
        for j in range(l):
            for k in range(n):
                C[i][j]+=A[i][k]*B[k][j]
                C[i][j]%=MOD
    return C

def resolve():
    n,k=map(int,input().split())
    A=[list(map(int,input().split())) for _ in range(n)]
    dp=[[1]]*n
    A=power(A,k,mat_mul)
    dp=mat_mul(A,dp)
    res=0
    for i in range(n): res+=dp[i][0]
    print(res%MOD)
resolve()

Solution (Python):
mod = 10**9 + 7

N, K = map(int,input().split())
A = [list(map(int,input().split())) for _ in range(N)]

def dot(X, Y):
    x = len(X)
    y = len(Y[0])
    z = len(Y)
    C = [[0] * y for _ in range(x)]
    for i in range(x):
        for j in range(y):
            for k in range(z):
                C[i][j] += X[i][k] * Y[k][j]
                C[i][j] %= mod
    return C

B = [[0] * N for _ in range(N)]
for k in range(N):
    B[k][k] = 1

while K > 0:
    if K & 1:
        B = dot(B, A)
    A = dot(A, A)
    K >>= 1

ans = 0
for i in range(N):
    for j in range(N):
        ans += B[i][j]
        ans %= mod

print(ans)

Solution (Python):
mod=1000000007


def dot(a,b):
    ret=[[0]*n for i in range(n)]
    for i in range(n):
        for j in range(n):
            ret[i][j]=sum(a[i][k]*b[k][j] for k in range(n))%mod
    return ret


n,k=map(int,input().split())
dp=[list(map(int,input().split())) for i in range(n)]
ans=[[0]*n for i in range(n)]
for i in range(n):
    ans[i][i]=1

while k:
    k,i=divmod(k,2)
    if i:
        ans=dot(ans,dp)
    dp=dot(dp,dp)
print(sum(sum(a)%mod for a in ans)%mod)

Solution (Python):
import math
MOD = int(1e9+7)


def main():
    N, K = map(int, input().split())
    A = [list(map(int, input().split())) for _ in range(N)]
    x = int(math.log(K, 2))
    dp = [[[0] * N for _ in range(N)] for _ in range(x+1)]
    dp[0] = A
    for d in range(x):
        for i in range(N):
            for j in range(N):
                for k in range(N):
                    dp[d+1][i][j] = (dp[d+1][i][j] +
                                     dp[d][i][k] * dp[d][k][j]) % MOD
    ans = [1] * N
    for d in range(x+1):
        if K & (2 ** d) == 0:
            continue
        prev = list(ans)
        ans = [0] * N
        for i in range(N):
            for j in range(N):
                ans[j] = (ans[j] + dp[d][i][j] * prev[i]) % MOD

    print(sum(ans) % MOD)


if __name__ == "__main__":
    main()

Solution (Python):
import sys
sys.setrecursionlimit(10000000)
def input():
    return sys.stdin.readline()[:-1]
from collections import *
import math
from bisect import *
import functools
INF = float('inf')
mod = 10**9+7

def mat(A, B):
    n = len(A)
    C = [[0] * n for i in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                C[i][j] += A[i][k]*B[k][j]
            C[i][j] %= mod
    return C

N, K = map(int, input().split())
a = [[int(x) for x in input().split()] for i in range(N)]
b = [[0] * N for i in range(N)]
for i in range(N):
    b[i][i] = 1
while K > 0:
    if K % 2:
        b = mat(b, a)
    a = mat(a, a)
    K //= 2
print(sum([sum(_b) for _b in b])%mod)

Solution (Python):


MOD=10**9+7

def matmul(A,B,M): 
    I,J,K = len(A), len(B[0]), len(B)  
    res = [[0] * J for _ in range(I)]
    for i in range(I):
        for k in range(K):
            for j in range(J):
                res[i][j] += A[i][k]*B[k][j]
                res[i][j] %= M
    return res

def matpow(A,p,M): 
    if p % 2:
        return matmul(A, matpow(A,p-1,M), M)
    elif p > 0:
        b = matpow(A, p//2, M)
        return matmul(b, b, M)
    else:
        return [[1 if i == j else 0 for j in range(len(A))] for i in range(len(A))]






n,k = [int(i) for i in input().split()]
A = [[int(i) for i in input().split()] for _ in range(n)]

ans = matpow(A,k,MOD)
print(sum(sum(r) % MOD for r in ans) % MOD)

Solution (Python):


import sys

class RWalk:
    def solve(self):

        MOD = 10 ** 9 + 7

        def mult(A, B):
            n, m, p = len(A), len(A[0]), len(B[0])
            assert(len(B) == m)
            C = [[0] * p for _ in range(n)]
            for i in range(n):
                for k in range(m):
                    a = A[i][k]
                    row = B[k]
                    for j in range(p):
                        C[i][j] = (C[i][j] + a * row[j]) % MOD

            return C

        def pow(A, k):
            n, m = len(A), len(A[0])
            if k == 0:
                return [[1 if i == j else 0 for j in range(m)] for i in range(n)]

            if k % 2 == 0:
                return pow(mult(A, A), k // 2)

            return mult(A, pow(A, k - 1))


        n, k = [int(_) for _ in input().split()]
        a = []
        for i in range(n):
            a.append([int(_) for _ in input().split()])

        print(sum([sum(line) for line in pow(a, k)]) % MOD)

solver = RWalk()
input = sys.stdin.readline

solver.solve()

Solution (Python):


import sys

class RWalk:
    def solve(self):

        MOD = 10 ** 9 + 7

        def mult(A, B):
            n, m, p = len(A), len(A[0]), len(B[0])
            assert(len(B) == m)
            C = [[0] * p for _ in range(n)]
            for i in range(n):
                for j in range(p):
                    for k in range(m):
                        C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD

            return C

        def pow(A, k):
            n, m = len(A), len(A[0])
            if k == 0:
                return [[1 if i == j else 0 for j in range(m)] for i in range(n)]

            if k % 2 == 0:
                return pow(mult(A, A), k // 2)

            return mult(A, pow(A, k - 1))


        n, k = [int(_) for _ in input().split()]
        a = []
        for i in range(n):
            a.append([int(_) for _ in input().split()])

        print(sum([sum(line) for line in pow(a, k)]) % MOD)

solver = RWalk()
input = sys.stdin.readline

solver.solve()

Solution (Python):
MOD = 10**9+7

def mat_dot(one, two, mod):
	return [[sum([(one[i][k]*two[k][j])%mod for k in range(max(len(two), len(one[0])))])%mod for j in range(len(two[0]))] for i in range(len(one))]

def mat_pow(mat, exp, mod):
	size = len(mat)
	res = [[0 for _ in range(size)] for _ in range(size)]
	for i in range(size):
		res[i][i] = 1
	cnt = 0
	while (1<<cnt) <= exp:
		if (exp>>cnt)&1:
			res = mat_dot(res, mat, mod)
		mat = mat_dot(mat, mat, mod)
		cnt += 1
	return res

n, k = map(int, input().split())
a = [list(map(int, input().split())) for _ in range(n)]
b = mat_pow(a, k, MOD)
print(sum([sum(x) for x in b]) % MOD)

Solution (Python):

n, k = [ int(v) for v in input().split() ]
path_mat = [ [ int(v) for v in input().split() ] for i in range(n) ]

mod = 10**9+7

def mat_mul(inmat1,inmat2,m):
    outmat = [ [ 0 for j in range(m) ] for i in range(m) ]

    for i in range(m):
        for j in range(m):
            for g in range(m):
                outmat[i][j] += inmat1[i][g] * inmat2[g][j]
            outmat[i][j] %= mod
    return outmat


ans_mat = [ [ 1 if i == j else 0 for j in range(n) ] for i in range(n) ]

while k != 0:
    if k % 2 == 1:
        ans_mat = mat_mul(ans_mat,path_mat,n)
    k = k // 2
    path_mat = mat_mul(path_mat,path_mat,n)
    

print(sum([ sum(i) for i in ans_mat ])%mod)

Solution (Python):
import numpy as np

N, K = list(map(int, input().split()))
a = [list(map(int, input().split())) for _ in range(N)]
MOD = 10 ** 9 + 7

A = np.array(a, dtype = 'O')
A = A.T

ans = np.eye(N, dtype = 'O')
while K:
  if K & 1:
    ans = np.dot(ans, A)
    ans %= MOD
  A = np.dot(A, A)
  A %= MOD
  K >>= 1

print(np.sum(abs(ans)) % MOD)

Solution (Python):
import sys
import numpy as np
input = sys.stdin.readline
sys.setrecursionlimit(10 ** 7)

MOD = 10 ** 9 + 7
N,K = map(int,input().split())
A = np.array([input().split() for _ in range(N)],dtype=np.int64)

def mult1(A,B):
    C = np.zeros_like(A)
    for i in range(N):
        for j in range(N):
            C[i,j] = (A[i,:] * B[:,j] % MOD).sum() % MOD
    return C

def mult(A,B):
    
    M = 1 << 16
    A1,A2 = A//M,A%M
    B1,B2 = B//M,B%M
    X,Y,Z = (np.dot(P,Q) % MOD for P,Q in ((A1,B1),(A1-A2,B1-B2),(A2,B2)))
    return ((X << 32) + (X+Z-Y << 16) + Z) % MOD

def power(A,k):
    if k == 1:
        return A
    X = power(A,k//2)
    X = mult(X,X)
    return mult(A,X) if k&1 else X

answer = power(A,K).sum() % MOD
print(answer)

Solution (Python):
N,K=list(map(int,input().strip().split(' ')))
A=[]
for _ in range(N):
    A+=[list(map(int,input().strip().split(' ')))]
dic={}
MOD=10**9+7
power=1
while 2**power<=K:
    power+=1

dp=[[[0 for j in range(N)] for i in range(N)] for p in range(power+1) ]
for p in range(0,power+1):
    
    if p==0:
        for i in range(N):
            for j in range(N):
                if A[i][j]==1:
                    dp[p][i][j]=1
    else:
        for i in range(N):
            for j in range(N):
                for k in range(N):
                    dp[p][i][j]+=(dp[p-1][i][k]*dp[p-1][k][j])%MOD
                    dp[p][i][j]%=MOD

    
rep=bin(K)[2:][::-1]
final_ans=[[0 for i in range(N)] for j in range(N)]
need=[]
for i in range(len(rep)):
    if rep[i]=="1":
        need+=[i]

final=[[dp[need[0]][i][j] for j in range(N)] for i in range(N)]
for i in range(1,len(need)):
    
    temp_ans=[[0 for i in range(N)] for j in range(N)]
    for x in range(N):
        for y in range(N):
            temp=0
            for z in range(N):
                temp+=final[x][z]*dp[need[i]][z][y]
                temp%=MOD
            temp_ans[x][y]=temp
    final=[[temp_ans[i][j] for j in range(N)] for i in range(N)]
            
ans=0
for i in range(N):
    for j in range(N):
        ans+=final[i][j]
        ans%=MOD
print(ans)

Solution (Python):
from collections import defaultdict,deque
import sys,heapq,bisect,math,itertools,string,queue,copy,time
sys.setrecursionlimit(10**8)
INF = float('inf')
mod = 10**9+7
eps = 10**-7
def inp(): return int(input())
def inpl(): return list(map(int, input().split()))
def inpl_str(): return list(input().split())

N,K = inpl()
aa = [inpl() for _ in range(N)]

def Matdot(A,B):
    C = [[0]*N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            for k in range(N):
                C[i][j] += A[i][k]*B[k][j]
            C[i][j] %= mod
    return C

def power(aa,k):
    if k == 1:
        return aa
    elif k%2 == 0:
        tmp = power(aa,k//2)
        return Matdot(tmp,tmp)
    else:
        tmp = power(aa,(k-1)//2)
        return Matdot(Matdot(tmp,tmp),aa)


aa =  power(aa,K)
ans = 0
for a in aa:
    ans += sum(a)
    ans %= mod
print(ans)

Solution (Python):
import sys
import numpy as np
input = sys.stdin.readline

MOD=10**9+7
N, K = map(int, input().split())
grid = np.empty((0,N), dtype='O')
for i in range(N):
  array = list(map(int, input().strip().split()))
  array2 = np.array(array)
  grid = np.append(grid,[array2],axis=0)

dp=[[grid]]
gridsave=grid
i=1
while 2**i<=K:
  
  
  gridsave=np.dot(dp[i-1][0],dp[i-1][0])
  gridsave=np.mod(gridsave,MOD)
  
  dp.append([gridsave])
  i+=1
i-=1
K-=2**i
i-=1

grid1=dp[len(dp)-1][0]


while K>0:
  if K>=2**i:
    
    grid1=np.dot(grid1,dp[i][0])
    grid1=np.mod(grid1,MOD)
    K-=2**i
    
  i-=1
  
ans=0

for i in range(N):
  for j in range(N):
    ans=(ans+grid1[i][j])%MOD
print(ans)

Solution (Python):
import sys
import numpy as np
input = sys.stdin.readline

MOD=10**9+7
N, K = map(int, input().split())
grid = np.empty((0,N), dtype='O')
for i in range(N):
  array = list(map(int, input().strip().split()))
  array2 = np.array(array)
  grid = np.append(grid,[array2],axis=0)

dp=[[grid]]
gridsave=np.eye(N,dtype="O")

i=1
while K:
  if K & 1:
    gridsave=np.dot(grid,gridsave) % MOD
  grid=np.dot(grid,grid) % MOD
  K>>=1
ans=np.sum(gridsave)%MOD

print(ans)

Solution (Python):
def mat_dot(A, B, mod):
    assert len(A[0]) == len(B), 'invalid_size'
    
    L = len(A)
    M = len(A[0])
    N = len(B[0])
    
    res = [[0]*N for _ in range(L)]
    
    for i in range(L):
        for j in range(N):
            a = 0
            for k in range(M):
                a = (a+A[i][k]*B[k][j]) % mod
            res[i][j] = a
            
    return res

def mat_pow(A, x, mod):
    N = len(A)
    res = [[0]*N for _ in range(N)]
    for i in range(N):
        res[i][i] = 1
    for i in range(x.bit_length()):
        if 2**i & x:
            res = mat_dot(res, A, mod)
        A = mat_dot(A, A, mod)
    return res

N, K = map(int, input().split())

A = [list(map(int, input().split())) for _ in range(N)]
A = list(map(list, zip(*A)))
mod = 10**9+7
B = mat_pow(A, K, mod)
ans = 0
for i in range(N):
    for j in range(N):
        ans = (ans + B[i][j]) % mod
print(ans)

Solution (Python):


MOD=10**9+7

def matmul(A, B, mod): 
    res = [[0] * len(B[0]) for _ in range(len(A))]
    for i, resi in enumerate(res):
        for k, aik in enumerate(A[i]):
            for j,bkj in enumerate(B[k]):
                resi[j] += aik*bkj
                resi[j] %= mod
    return res

def matpow(A, p, M): 
    if p % 2 == 1:
        return matmul(A, matpow(A, p - 1, M), M)
    elif p > 0:
        b = matpow(A, p // 2, M)
        return matmul(b, b, M)
    else:
        return [[1 if i == j else 0 for j in range(len(A))] for i in range(len(A))]


n,k = [int(i) for i in input().split()]
A = [[int(i) for i in input().split()] for _ in range(n)]

ans = matpow(A,k,MOD)
print(sum(sum(r) % MOD for r in ans) % MOD)

Solution (Python):


MOD=10**9+7

def dot(a, b):
    res = [[0] * n for _ in range(n)]
    for i in range(n):
        for k, aik in enumerate(a[i]):
            for j,bkj in enumerate(b[k]):
                res[i][j] += aik*bkj
                res[i][j] %= MOD
    return res


n,k = [int(i) for i in input().split()]
A = [[int(i) for i in input().split()] for _ in range(n)]

ans = [[0] * n for _ in range(n)]
for i in range(n):
    ans[i][i] = 1

while k:
    k, i = divmod(k, 2)
    if i:
        ans = dot(ans, A)
    A = dot(A, A)

res = 0
for i in range(n):
  for j in range(n):
    res += ans[i][j]
    res %= MOD
print(res)

Solution (Python):


MOD=10**9+7

def dot(a, b):
    res = [[0] * n for _ in range(n)]
    for i in range(n):
        for k, aik in enumerate(a[i]):
            for j,bkj in enumerate(b[k]):
                res[i][j] += aik*bkj
                res[i][j] %= MOD
 
    return res


n,k = [int(i) for i in input().split()]
A = [[int(i) for i in input().split()] for _ in range(n)]

ans = [[0] * n for _ in range(n)]
for i in range(n):
    ans[i][i] = 1

while k:
    k, i = divmod(k, 2)
    if i:
        ans = dot(ans, A)
    A = dot(A, A)

print(sum(sum(A) % MOD for A in ans) % MOD)

Solution (Python):
n, k = map(int, input().split())
dp = [list(map(int, input().split())) for i in range(n)]
MOD = 10**9 + 7

def _mul(A, B, MOD):
    C = [[0] * len(B[0]) for i in range(len(A))]
    for i in range(len(A)):
        for k in range(len(B)):
            for j in range(len(B[0])):
                C[i][j] = (C[i][j] + A[i][k]*B[k][j]) % MOD
    return C


def pow(A, n, MOD):
    B = [[0] * len(A) for i in range(len(A))]
    for i in range(len(A)):
        B[i][i] = 1
    while n > 0:
        if n & 1:
            B = _mul(A, B, MOD)
        A = _mul(A, A, MOD)
        n  = n // 2
    return B

ans_matrix = pow(dp, k, MOD)
ans = 0
for i in range(n):
    ans += sum(ans_matrix[i])
    ans %= MOD
print(ans)

Solution (Python):
import numpy as np
import sys

input = sys.stdin.readline

md = 10 ** 9 + 7
n, k = map(int, input().split())
a = np.array([list(map(int, input().split())) for _ in range(n)], dtype="O")
c = np.eye(n, dtype="O")

while k:
    if k & 1:
        c = np.dot(a, c) % md
    a = np.dot(a, a) % md
    k >>= 1
print(np.sum(c) % md)

Solution (Python):
import sys

input = sys.stdin.readline


def dot(a, b):
    c = [[0] * n for _ in range(n)]
    for i, ai in enumerate(a):
        for j, bj in enumerate(zip(*b)):
            c[i][j] = sum(ak * bk for ak, bk in zip(ai, bj)) % md
    return c


def pwr(a, k):
    c = None
    while c == None:
        if k & 1:
            c = a
        a = dot(a, a)
        k >>= 1

    while k:
        if k & 1:
            c = dot(a, c)
        a = dot(a, a)
        k >>= 1
    return c


md = 10 ** 9 + 7
n, k = map(int, input().split())
a = [list(map(int, input().split())) for _ in range(n)]
print(sum([sum(pk) % md for pk in pwr(a, k)]) % md)