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):
from numpy import array, eye

mod = 10 ** 9 + 7

N, K, *A = map(int, open(0).read().split())
A = array(A, dtype=object).reshape(N, N)

ans = eye(N, dtype=object)
while K:
    if K & 1:
        ans = (ans.dot(A)) % mod
    A = A.dot(A) % mod
    K //= 2

print(ans.sum() % mod)

Solution (Python):
n, kk = map(int, input().split())
AA = [list(map(int, input().split())) for _ in range(n)]
mod = 10 ** 9 + 7




def matrix_product(A, B):
    res = [[0] * n for _ in range(n)]
    for i in range(n):
        for k in range(n):
            for j in range(n):
                res[i][j] += A[i][k] * B[k][j]
                res[i][j] %= mod
    return res


def matrix_exponents(A, t):
    res = [[0] * n for _ in range(n)]
    for i in range(n):
        res[i][i] = 1

    while t:
        if t & 1 == 1:
            res = matrix_product(res, A)
        A = matrix_product(A, A)
        t //= 2
    return res


temp = matrix_exponents(AA, kk)


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

Solution (Python):
mod = 10**9+7
n, K = map(int, input().split())
a = [[int(x) for x in input().split()] for _ in range(n)]

A = [a]

i = 1
while 2**i <= K:
    mtrix = []
    for j in range(n):
        raw = []
        for l in range(n):
            tmp = 0
            for k in range(n):
                tmp += A[-1][j][k] * A[-1][k][l]
            raw.append(tmp % mod)
        mtrix.append(raw)
    A.append(mtrix)
    i += 1

bit_max = i

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

for i in range(bit_max):
    if (K >> i) & 1:
        mtrix = []
        for j in range(n):
            raw = []
            for l in range(n):
                tmp = 0
                for k in range(n):
                    tmp += ans[j][k] * A[i][k][l]
                raw.append(tmp % mod)
            mtrix.append(raw)
        ans = mtrix

ans = sum(map(lambda x: sum(x) % mod, ans))

print(ans % mod)

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

def mat_pow(A, n):
    res = [[0]*len(A) for _ in range(len(A))]
    
    for i in range(len(A)):
        res[i][i] = 1
    
    while n>0:
        if n&1:
            res = mat_mul(res, A)
        
        n >>= 1
        A = mat_mul(A, A)
    
    return res

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

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

print(ans)

Solution (Python):
import sys
def input(): return sys.stdin.readline().rstrip()
def multi_mat(M1,M2,n,mod=10**9+7):
    ans=[[0]*n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                ans[i][j]=(ans[i][j]+M1[i][k]*M2[k][j])%mod
    return ans
def pow_mat(M,k,n,mod=10**9+7):
    ans=[[0]*n for _ in range(n)]
    for i in range(n):ans[i][i]=1
    while k:
        if k%2:ans=multi_mat(ans,M,n,mod)
        M=multi_mat(M,M,n,mod)
        k//=2
    return ans

def main():
    n,k=map(int,input().split())
    A=[[int(_) for _ in input().split()] for i in range(n)]
    mod=10**9+7
    Ax=pow_mat(A,k,n,mod)
    print(sum(sum(x)%mod for x in Ax)%mod)

if __name__=='__main__':
    main()

Solution (Python):

def matmul(mt1, mt2, n, modp):
    ans = [[0 for _ in range(n)] for _ in range(n)]

    for i in range(n):
        for j in range(n):
            for k in range(n):
                ans[i][j] += mt1[i][k] * mt2[k][j]
                ans[i][j] %= modp
    return ans


def pow_matrix(mt, k, n, modp):
    ans = [[0 for _ in range(n)] for _ in range(n)]
    for i in range(n):
        ans[i][i] = 1

    while k:
        if k % 2:
            ans = matmul(ans, mt, n, modp)
        mt = matmul(mt, mt, n, modp)
        k //= 2
    return ans


def submit():
    n, k = map(int, input().split())
    a = [list(map(int, input().split())) for _ in range(n)]
    modp = 10 ** 9 + 7
    ak = pow_matrix(a, k, n, modp)
    
    ans = 0
    for ar in ak:
        for ac in ar:
            ans += ac
            ans %= modp
    print(ans)

submit()

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF = float("inf")
MOD = 10**9 + 7 
input = lambda:sys.stdin.readline().rstrip()

def resolve():
    n, p = map(int, input().split())
    
    base = [list(map(int, input().split())) for _ in range(n)]
    dp = [[(i == j) + 0 for j in range(n)] for i in range(n)]
    while p:
        if p & 1:
            dp = [[sum(dp[i][k] * base[k][j] % MOD for k in range(n)) % MOD for j in range(n)] for i in range(n)]
        base = [[sum(base[i][k] * base[k][j] % MOD for k in range(n)) % MOD for j in range(n)] for i in range(n)]
        p >>= 1

    print(sum(dp[i][j] for i in range(n) for j in range(n)) % MOD)
resolve()

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF = float("inf")
MOD = 10**9 + 7 
input = lambda:sys.stdin.readline().rstrip()

def resolve():
    n, p = map(int, input().split())
    
    base = [list(map(int, input().split())) for _ in range(n)]
    dp = [[(i == j) + 0 for j in range(n)] for i in range(n)]
    while p:
        if p & 1:
            dp = [[sum(dp[i][k] * base[k][j] for k in range(n)) % MOD for j in range(n)] for i in range(n)]
        base = [[sum(base[i][k] * base[k][j] for k in range(n)) % MOD for j in range(n)] for i in range(n)]
        p >>= 1

    print(sum(dp[i][j] for i in range(n) for j in range(n)) % MOD)
resolve()

Solution (Python):
import sys
sys.setrecursionlimit(100000000)
input = sys.stdin.readline
MOD = 10 ** 9 + 7

class Matrix(list):

    def mul(self,B,mod = 0):
        C = Matrix([[0] * len(B[0]) for _ in range(len(self))])
        if mod == 0:
            for i in range(len(self)):
                for k in range(len(B)):
                    for j in range(len(B[0])):
                        C[i][j] += self[i][k] * B[k][j]
            
            return C
        else:
            for i in range(len(self)):
                for k in range(len(B)):
                    for j in range(len(B[0])):
                        C[i][j] = (C[i][j] + self[i][k] * B[k][j]) % mod
            
            return C

    
    def pow(self,n,mod = 0):
        if len(self) != len(self[0]):
            raise ValueError('the size of matrix is different')

        if mod == 0:
            B = Matrix([[0] * len(self) for _ in range(len(self))])
            A = self
            for i in range(len(self)):
                B[i][i] = 1
            
            while n > 0:
                if (n&1):
                    B = B.mul(A)
                A = A.mul(A)
                n //= 2
            
            return B

        else:
            B = Matrix([[0] * len(self) for _ in range(len(self))])
            A = self
            for i in range(len(self)):
                B[i][i] = 1
            
            while n > 0:
                if (n&1):
                    B = B.mul(A,mod)
                A = A.mul(A,mod)
                n //= 2
            
            return B

def main():
    N,K = map(int,input().split())
    A = Matrix([list(map(int,input().split())) for _ in range(N)])
    B = A.pow(K,MOD)
    ans = sum(sum(B[i])%MOD for i in range(N))%MOD
    print(ans)
if __name__ == '__main__':
    main()

Solution (Python):
import sys
input = sys.stdin.readline
mod = 10**9 + 7
n,k = map(int,input().split())
x = [list(map(int,input().split())) for i in range(n)]
ans = [[0]*i +[1]+[0]*(n-i-1) for i in range(n)]

def dot(A,B):
    res = [[0]*n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                res[i][j] += A[i][k]*B[k][j]
                res[i][j] %= mod
    return res

while k > 0:
    if k&1:
        ans = dot(ans, x)
    k = k >> 1
    x = dot(x,x)
res = 0
for i in range(n):
    for j in range(n):
        res += ans[i][j]
        res %= mod
print(res)

Solution (Python):
import numpy as np
mod = 10 ** 9 + 7


def matpow(A, x):
    if x == 1:
        return A % mod
    elif x % 2 == 0:
        return matpow(np.linalg.matrix_power(A, 2) % mod, x // 2)
    else:
        ans = matpow(np.linalg.matrix_power(A, 2) % mod, x // 2)
        ans = ans.dot(A) % mod
        return ans


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

A = np.array(A, dtype=object)
print(np.sum(matpow(A, K) % mod) % mod)

Solution (Python):




P=10**9+7
N,K=map(int,input().split())
G=[]
for index in range(N):
    G.append(list(map(int,input().split())))






def mul(x,y,p): 
    res=[[0]*(len(y[0])) for _ in range(len(x))]
    
    for index1 in range(len(x)):
        for index2 in range(len(y[0])):
            for index3 in range(len(y)):  
                res[index1][index2]=(res[index1][index2]+x[index1][index3]*y[index3][index2])%p
    return res


def pow(x,y,p):
    res=[[0]*(len(x)) for _ in range(len(x))]
    for index1 in range(len(x)):
        res[index1][index1]=1
    
    while y>0:
        if y&1==1:
            res=mul(res,x,p)
        x=mul(x,x,p)
        y=y>>1
    
    return res




DP_K=pow(G,K,P)

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

print(ans)

Solution (Python):
def product(A, B):
    C = [[0]*n for i in range(n)]
    for i in range(n):
        for j in range(n):
            C[i][j] = sum([A[i][k]*B[k][j] for k in range(n)]) % p
    return C

def matpow(A, k):
    if k == 1:
        return A
    if k % 2 == 0:
        B = matpow(A, k//2)
        return product(B, B)
    return product(A, matpow(A, k-1))

n, k = map(int, input().split())
p = 10**9+7
A = [list(map(int, input().split())) for i in range(n)]
B = matpow(A, k)
print(sum([sum(B[i]) for i in range(n)]) % p)

Solution (Python):
import sys

sys.setrecursionlimit(10 ** 7)
rl = sys.stdin.readline


def matmul(a, b, mod):
    N0 = len(a)
    N1 = len(b[0])
    N2 = len(a[0])
    res = [[0] * N1 for _ in range(N0)]
    for i in range(N0):
        for j in range(N1):
            tmp = 0
            for k in range(N2):
                tmp = (tmp + a[i][k] * b[k][j]) % mod
            res[i][j] = tmp
    return res


def matpow(a, k, mod):
    N = len(a)
    res = [[0] * N for _ in range(N)]
    for i in range(N):
        res[i][i] = 1
    while k:
        if k & 1:
            res = matmul(res, a, mod)
        a = matmul(a, a, mod)
        k >>= 1
    return res


def solve():
    N, K = map(int, rl().split())
    a = [list(map(int, rl().split())) for _ in range(N)]
    MOD = 10 ** 9 + 7
    
    a = list(map(list, zip(*a)))
    mat = matpow(a, K, MOD)
    ans = 0
    for i in range(N):
        for j in range(N):
            ans = (ans + mat[i][j]) % MOD
    print(ans)


if __name__ == '__main__':
    solve()

Solution (Python):
N, K = map(int, input().split())
A = [list(map(int, input().split())) for i in range(N)]
mod = 10**9+7

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

now = [[1 if i == j else 0 for i in range(N)] for j in range(N)]        

while K:
    if K&1 == 1:
        now = matrixmul(now, A)
    K //= 2
    A = matrixmul(A, A)

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

Solution (Python):
import numpy as np
mod=10**9+7
def solve():
    n,k = map(int,input().split())
    a = np.array([list(map(int,input().split())) for i in range(n)],dtype=np.int64)
    print(modpow(a,k).sum()%mod)
def mult(A,B):
    N=len(A)
    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 modpow(A,k):
    n=len(A)
    res=np.eye(n,dtype=np.int64)
    for i in range(k.bit_length()):
        if k>>i&1:
            res=mult(A,res)
        A=mult(A,A)
    return res
solve()

Solution (Python):


import sys
input = sys.stdin.readline

def dot(A,B,MOD=1000000007):
    N,M,L = len(A),len(A[0]),len(B[0])
    res = [[0]*L for i in range(N)]
    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    return res

def matPow(A,x,MOD=1000000007):
    N = len(A)
    res = [[0]*N for i in range(N)]
    for i in range(N):
        res[i][i] = 1
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = dot(res,A)
        A = dot(A,A)
    return res

def main():
    n,k = map(int,input().split())
    mod = 10**9+7
    a = [list(map(int,input().split())) for i in range(n)]
    dp = matPow(a,k)
    ans = 0
    for i in range(n):
        for j in range(n):
            ans += dp[i][j]
    print(ans%mod)

main()

Solution (Python):


n,k=map(int,input().split())
a=[list(map(int,input().split())) for i in range(n)]
dp=[[0]*n for i in range(n)]
mod=10**9+7
for i in range(n):
    for j in range(n):
        dp[i][i]=1
def dot(A,B):
    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
while k:
    if k&1:
        dp=dot(dp,a)
    a=dot(a,a)
    k>>=1
ans=0
for i in dp:
    ans+=sum(i)
print(ans%mod)

Solution (Python):

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))

def main():
    mod=10**9+7
    N,K=MI()
    a=[[]for _ in range(N)]
    for i in range(N):
        a[i]=LI()
        
    def matDot(A,B,MOD):
        N,M,L = len(A),len(A[0]),len(B[0])
        
        res = [[0]*L for i in range(N)]

        for i in range(N):
            for j in range(L):
                s = 0
                for k in range(M):
                    s = (s + A[i][k]*B[k][j]) % MOD
                res[i][j] = s
        
        return res

    def matPow(A,x,MOD):
        N = len(A)
        res = [[0]*N for i in range(N)]
        
        for i in range(N):
            res[i][i] = 1
        
        for i in range(x.bit_length()):
            if (x>>i) & 1:
                res = matDot(res,A,MOD)
            A = matDot(A,A,MOD)
    
        return res
    
    to=matPow(a,K,mod)
    ans=0
    for i in range(N):
        for j in range(N):
            ans=(ans+to[i][j])%mod
    print(ans)
    
        

main()

Solution (Python):

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))

def main():
    mod=10**9+7
    N,K=MI()
    a=[[]for _ in range(N)]
    for i in range(N):
        a[i]=LI()
    
    adj=[[]for _ in range(N)]
    for i in range(N):
        for j in range(N):
            if a[i][j]==1:
                adj[i].append(j)
    
        
    M=1
    mm=1
    for i in range(70):
        mm*=2
        M+=1
        if mm>=K:
            break
    
    dab=[[[0]*N for _ in range(N)]for i in range(M)]
    
    for j in range(N):
        for kk in range(len(adj[j])):
            k=adj[j][kk]
            dab[0][j][k]=1
    
    for i in range(M-1):
        for j in range(N):
            for k in range(N):
                cnt=dab[i][j][k]
                for l in range(N):
                    dab[i+1][j][l]+=cnt*dab[i][k][l]
                    dab[i+1][j][l]%=mod
                    
    ans=[1]*N
    for mm in range(M):
        if K>>mm&1:
            temp=[0]*N
            for i in range(N):
                cnt=ans[i]
                for j in range(N):
                    temp[j]+=dab[mm][i][j]*cnt
                    
            for i in range(N):
                ans[i]=temp[i]
                
    res=0
    for i in range(N):
        res+=ans[i]
        res%=mod
    print(res)
        

main()

Solution (Python):
def matmul(p,q):
  c = [[0 for i in range(n)] for j in range(n)]
  for i in range(n):
    for j in range(n):
      for l in range(n):
        c[i][j] += p[i][l]*q[l][j]
        c[i][j] %= mod
  return c
n,k = map(int,input().split())
a = [list(map(int,input().split())) for i in range(n)]
mod = 10**9+7
als = [a]
for i in range(100):
  alsnxt = matmul(als[-1],als[-1])
  als.append(alsnxt)
ansmat = [[0 for i in range(n)] for j in range(n)]
for i in range(n):
  ansmat[i][i] = 1
for i in range(100):
  if k&1<<i:
    ansmat = matmul(ansmat,als[i])
ans = 0
for i in range(n):
  for j in range(n):
    ans += ansmat[i][j]
    ans %= mod
print(ans)

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

mod=10**9+7

dp=[[[0 for i in range(0,N)] for j in range(0,N)] for k in range(0,62)]

for j in range(0,N):
    for k in range(0,N):
        dp[0][j][k]=A[j][k]

for i in range(1,62):
    for j in range(0,N):
        for k in range(0,N):
            dp[i][j][k]=sum(dp[i-1][j][l]*dp[i-1][l][k] for l in range(0,N))%mod

ans=[[[0 for i in range(0,N)] for j in range(0,N)] for k in range(0,63)]

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

for k in range(0,62):
    if K>>k & 1==1:
        for i in range(0,N):
            for j in range(0,N):
                ans[k+1][i][j]=sum(dp[k][i][l]*ans[k][l][j] for l in range(0,N))%mod

    else:
        for i in range(0,N):
            for j in range(0,N):
                ans[k+1][i][j]=ans[k][i][j]

print(sum(ans[62][i][j] for i in range(0,N) for j in range(0,N))%mod)

Solution (Python):
import numpy as np

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

a = np.matrix(a, dtype=object)
ans = np.matrix(np.diag([1]*N), dtype=object)

while K:
    K, m = divmod(K, 2)
    if m:
        ans = ans * a
        ans = np.mod(ans, MOD)
    a = a*a
    a = np.mod(a, MOD)

print(ans.sum() % MOD)

Solution (Python):
n, k = [int(i) for i in input().split()]
mod = 10 ** 9 + 7

a = [[int(i) for i in input().split()]for _ in range(n)]
b = [[0] * n for _ in range(n)]
for i in range(n):
    b[i][i] = 1


def mul(a, b):
    c = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            s = 0
            for k in range(n):
                s += a[i][k] * b[k][j]
                s %= mod
            c[i][j] = s
    return c


while k != 0:
    if k & 1 == 1:
        b = mul(b, a)
    a = mul(a, a)
    k //= 2

ans = sum(map(sum, b)) % mod

print(ans)

Solution (Python):
from numpy import*
(n,k),*a=[list(map(int,t.split()))for t in open(0)]
a=matrix(a,'O')
M=10**9+7
s=1
while k:s*=~k%2or a;k//=2;a=a*a%M
print(sum(s)%M)

Solution (Python):
M=10**9+7
(n,k),*a=[list(map(int,t.split()))for t in open(0)]
r=a if k&1else[[i==j for j in range(n)]for i in range(n)]
d=a
for i in range(1,k.bit_length()):
  p=[]
  for y in d:
    q=[0]*n
    for x,s in zip(y,d):
      for j,t in enumerate(s):
        q[j]+=t*x%M
    p+=q,
  d=p
  if 1<<i&k:
    p=[]
    for y in r:
      q=[0]*n
      for x,s in zip(y,d):
        for j,t in enumerate(s):
          q[j]+=t*x%M
      p+=q,
    r=p
print(sum(t for s in r for t in s)%M)

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

def main():
    N,K = map(int,input().split())
    MOD = 10**9+7
    
    a = [list(map(int,input().split())) for _ in range(N)]
    arr = np.array(a,dtype=object)

    use = np.eye(N,dtype=object)
    while K>0:
        if K&1:
            use = np.dot(use,arr)%MOD
        arr = np.dot(arr,arr)%MOD
        K //= 2

    print(np.sum(use)%MOD)

if __name__ == "__main__":
    main()

Solution (Python):
import sys, re
from collections import deque, defaultdict, Counter
from math import ceil, sqrt, hypot, factorial, pi, sin, cos, radians
from itertools import accumulate, permutations, combinations, product
from operator import itemgetter
from copy import deepcopy
from string import ascii_lowercase, ascii_uppercase, digits
from bisect import bisect, bisect_left
from fractions import gcd
from heapq import heappush, heappop
from functools import reduce, lru_cache
from pprint import pprint
import numpy as np
def input(): return sys.stdin.readline().strip()
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(): return list(map(int, input().split()))
def ZIP(n): return zip(*(MAP() for _ in range(n)))
sys.setrecursionlimit(10 ** 9)
INF = float('inf')
mod = 10 ** 9 + 7

def matpow(A, n):  
	B = np.identity(len(A), dtype=object)
	while n > 0:
		if n & 1:
			B = B.dot(A)
			B %= mod
		A = A.dot(A)
		A %= mod
		n >>= 1
	return B

N, K = MAP()
a = np.array([LIST() for _ in range(N)], dtype=object)

print(matpow(a, K).sum()%mod)

Solution (Python):
import sys
input = sys.stdin.readline
N, K = map(int, input().split())
mod = 10 ** 9 + 7
d = [list(map(int, input().split())) for _ in range(N)]
ln = K.bit_length()

def ArrayMultiple(a, b, mod):
  n = len(a)
  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

bn = list(bin(K)[2: ])
bn.reverse()
res = [[int(i == x) for i in range(N)] for x in range(N)]
for i in range(ln):
  if int(bn[i]):
    res = ArrayMultiple(res, d, mod)
  d = ArrayMultiple(d, d, mod)
rres = 0
for r in res:
  for x in r:
    rres += x
    rres %= mod
print(rres)

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

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

def power(base, exp, op=lambda x, y: x * y, ie=1):
    res = ie
    while exp:
        if exp & 1:
            res = op(res, base)
        base = op(base, base)
        exp >>= 1
    return res

matrix_id = [[i == j for j in range(N)] for i in range(N)]

def matrix_mul(A, B):
    res = [[0 for j in range(N)] for i in range(N)]
    for i in range(N):
        for j in range(N):
            for k in range(N):
                res[i][j] += A[i][k] * B[k][j]
                res[i][j] %= MOD
    return res

P = power(A, K, matrix_mul, matrix_id)

res = 0

for i in range(N):
    for j in range(N):
        res += P[i][j]
        res %= MOD

print(res)

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

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

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

def mod_pow(matrix, n, mod): 
    
    N = len(matrix[0])
    if n == 0:
        lst = [[1 if i == j else 0 for i in range(N)] for j in range(N)]
        return lst
    lst = [[0] * N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            for k in range(N):
                
                lst[i][j] += matrix[i][k] * matrix[k][j]
            lst[i][j] %= mod
    res = mod_pow(lst, n // 2, mod)

    lst = [[0] * N for _ in range(N)]
    if n % 2 == 1:
        for i in range(N):
            for j in range(N):
                for k in range(N):
                    lst[i][j] += res[i][k] * matrix[k][j]
                lst[i][j] %= mod
        return lst
    return res 

matrix = mod_pow(A, K, MOD)


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

ans %= MOD
print (ans)

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

def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

def main():
    n,k = map(int,input().split())
    mod = 10**9+7

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

    b = matPow(a,k,mod)

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

    print(ans)

if __name__ == "__main__":
    main()

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

def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

n,k = map(int,input().split())
mod = 10**9+7

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

b = matPow(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):
import sys
input = sys.stdin.readline

def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

n,k = map(int,input().split())
mod = 10**9+7

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

b = matPow(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):
import sys
input = sys.stdin.readline

def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

n,k = map(int,input().split())
mod = 10**9+7

a = [list(map(int,input().split())) for i in range(n)]

b = matPow(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):
def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (2**i) & x:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

n,k = map(int,input().split())
mod = 10**9+7

a = [list(map(int,input().split())) for i in range(n)]

b = matPow(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):
def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

n,k = map(int,input().split())
mod = 10**9+7

a = [list(map(int,input().split())) for i in range(n)]
a = list(map(list, zip(*a)))

b = matPow(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):
def matDot(A,B,MOD):
    N,M,L = len(A),len(A[0]),len(B[0])
    
    res = [[0]*L for i in range(N)]

    for i in range(N):
        for j in range(L):
            s = 0
            for k in range(M):
                s = (s + A[i][k]*B[k][j]) % MOD
            res[i][j] = s
    
    return res

def matPow(A,x,MOD):
    N = len(A)
    res = [[0]*N for i in range(N)]
    
    for i in range(N):
        res[i][i] = 1
    
    for i in range(x.bit_length()):
        if (x>>i) & 1:
            res = matDot(res,A,MOD)
        A = matDot(A,A,MOD)
    
    return res

n,k = map(int,input().split())
mod = 10**9+7

a = [list(map(int,input().split())) for i in range(n)]
b = matPow(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):
import numpy as np
n,k = map(int, input().split())
mod = 10**9 + 7
_a = [list(map(int, input().split())) for _ in range(n)]
a = np.array(_a, dtype = object)
e = np.eye(n, dtype = object)
while k:
  if k & 1:
    e = np.dot(e, a) % mod
  k >>= 1
  a = np.dot(a, a) % mod
print(e.sum() % mod)