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):
mod = 10**9+7
def mult(A, B):
    C = [[0 for i in range(len(A))] for i in range(len(A))]

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

def exp(A, p):
    ans = [[0 for i in range(len(A))] for i in range(len(A))]
    for i in range(len(A)):
        ans[i][i] = 1
    
    while p:
        if p%2 == 0:
            A = mult(A, A)
            p = p//2
        else:
            ans = mult(ans, A)
            p -= 1
    
    return ans

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

print(sum([sum(i)%mod for i in exp(grid, k)])%mod)

Solution (Python):
def mul(x,y):
    q=[[0]*n for i in xrange(n)]
    for i in xrange(n):
        for j in xrange(n):
            for k in xrange(n):
                q[i][j]+=x[i][k]*y[k][j]
    for i in xrange(n):
        for j in xrange(n):
            x[i][j]=q[i][j]%mod
n,k=map(int,raw_input().split())
a=[map(int,raw_input().split()) for _ in xrange(n)]
mod=pow(10,9)+7
z=[[int(i==j) for j in xrange(n)] for i in xrange(n)]
while k:
    if k%2:
        mul(z,a)
    mul(a,a)
    k/=2
print sum(sum(i) for i in z)%mod

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

mod = 10**9+7
n, k = map(int, input().split())
A = tuple(tuple(map(int, input().split())) for _ in range(n))
def matmul(X, Y):
  Z = [[0]*n for _ in range(n)]
  for i in range(n):
    for j in range(n):
      for k in range(n):
        Z[i][j] += X[i][k] * Y[k][j]
        Z[i][j] %= mod
  return Z
  
def matpow(k):
  if k == 1:
    return A
  if k%2:
    return matmul(A, matpow(k-1))
  X = matpow(k//2)
  return matmul(X, X)

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

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W)->np.ndarray:
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


N, K = ints()
A = read_matrix(N, N)

















def matmul_mod(X: np.ndarray, Y: np.ndarray, MOD: int=10**9 + 7)->np.ndarray:
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    YT = Y.T
    for i in range(X.shape[0]):  
        tmp = (YT * X[i, :]) % MOD
        ret[i, :] = tmp.sum(1) % MOD
    return ret










A_double = np.asanyarray(A, dtype=np.int64)
ans = np.eye(N, dtype=np.int64)
for j in range(int(K).bit_length()):
    if (K >> j) & 1:
        ans = matmul_mod(ans, A_double)
    A_double = matmul_mod(A_double, A_double)

print(ans.sum() % MOD)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W)->np.ndarray:
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


N, K = ints()
A = read_matrix(N, N)






def matmul_mod(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    for i in range(X.shape[0]):  
        for j in range(Y.shape[1]):  
            tmp = (X[i, :] * Y[:, j]) % MOD
            ret[i, j] = tmp.sum() % MOD
    return ret


def matmul_mod_try(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    YT = Y.T
    for i in range(X.shape[0]):  
        Xi = X[i, :]
        tmp = (YT * Xi) % MOD
        ret[i, :] = tmp.sum(1) % MOD
    return ret


def matmul_mod_fast(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    X3d = X[:, :, np.newaxis].transpose(0, 2, 1)
    Y3d = Y[:, :, np.newaxis].transpose(2, 1, 0)
    return ((X3d * Y3d) % MOD).sum(2) % MOD


A_double = np.asanyarray(A, dtype=np.int64)
ans = np.eye(N, dtype=np.int64)
for j in range(int(K).bit_length()):
    if (K >> j) & 1:
        ans = matmul_mod_fast(ans, A_double)
    A_double = matmul_mod_fast(A_double, A_double)

print(ans.sum() % MOD)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W)->np.ndarray:
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


N, K = ints()
A = read_matrix(N, N)






def matmul_mod(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    for i in range(X.shape[0]):  
        for j in range(Y.shape[1]):  
            tmp = (X[i, :] * Y[:, j]) % MOD
            ret[i, j] = tmp.sum() % MOD
    return ret


def matmul_mod_try(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    YT = Y.T
    for i in range(X.shape[0]):  
        Xi = X[i, :]
        tmp = (YT * Xi) % MOD
        ret[i, :] = tmp.sum(1) % MOD
    return ret


def matmul_mod_fast(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    X3d = X[:, :, np.newaxis].transpose(0, 2, 1)
    Y3d = Y[:, :, np.newaxis].transpose(2, 1, 0)
    return ((X3d * Y3d) % MOD).sum(2) % MOD


A_double = np.asanyarray(A, dtype=np.int64)
ans = np.eye(N, dtype=np.int64)
for j in range(int(K).bit_length()):
    if (K >> j) & 1:
        
        
        ans = matmul_mod_fast(ans, A_double)
    A_double = matmul_mod_fast(A_double, A_double)

print(ans.sum() % MOD)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W)->np.ndarray:
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


N, K = ints()
A = read_matrix(N, N)






def matmul_mod(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    for i in range(X.shape[0]):  
        for j in range(Y.shape[1]):  
            tmp = (X[i, :] * Y[:, j]) % MOD
            ret[i, j] = tmp.sum() % MOD
    return ret


def matmul_mod_try(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    YT = Y.T
    for i in range(X.shape[0]):  
        Xi = X[i, :]
        tmp = (YT * Xi) % MOD
        ret[i, :] = tmp.sum(1) % MOD
    return ret


A_double = np.asanyarray(A, dtype=np.int64)
ans = np.eye(N, dtype=np.int64)
for j in range(int(K).bit_length()):
    if (K >> j) & 1:
        
        
        ans = matmul_mod_try(ans, A_double)
    A_double = matmul_mod_try(A_double, A_double)

print(ans.sum() % MOD)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W)->np.ndarray:
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


N, K = ints()
A = read_matrix(N, N)






def matmul_mod(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    for i in range(X.shape[0]):  
        for j in range(Y.shape[1]):  
            tmp = (X[i, :] * Y[:, j]) % MOD
            ret[i, j] = tmp.sum() % MOD
    return ret


A_double = np.asanyarray(A, dtype=np.int64)
ans = np.eye(N, dtype=np.int64)
for j in range(int(K).bit_length()):
    if (K >> j) & 1:
        ans = matmul_mod(ans, A_double)
    A_double = matmul_mod(A_double, A_double)

print(ans.sum() % MOD)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)


def ints(): return np.fromstring(readline(), sep=' ', dtype=np.int64)


def read_matrix(H, W)->np.ndarray:
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


N, K = ints()
A = read_matrix(N, N)


from numba import njit


@njit('i8[:,:](i8[:,:],i8[:,:])', cache=True)
def matmul_mod(X: np.ndarray, Y: np.ndarray)->np.ndarray:
    MOD = 10**9 + 7
    
    ret = np.empty((X.shape[0], Y.shape[1]), dtype=np.int64)
    for i in range(X.shape[0]):  
        for j in range(Y.shape[1]):  
            tmp = (X[i, :] * Y[:, j]) % MOD
            ret[i, j] = tmp.sum() % MOD
    return ret


A_double = np.asanyarray(A, dtype=np.int64)
ans = np.eye(N, dtype=np.int64)
for j in range(int(K).bit_length()):
    if (K >> j) & 1:
        ans = matmul_mod(ans, A_double)
    A_double = matmul_mod(A_double, A_double)

print(ans.sum() % MOD)

Solution (Python):
import numpy as np
N,K=map(int,input().split())
a=[[int(i) for i in input().split()] for _ in range(N)]
a=np.array(a,dtype=np.int64)
MOD=10**9+7

def matmul(A,B,MOD):
    N=len(A)
    C=np.zeros((N,N),dtype=np.int64)

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

def mat_power(A,k,MOD):
    N=len(A)
    D=np.identity(N,dtype=np.int64)
    while k:
        if k & 1:
            D=matmul(D,A,MOD)
        A = matmul(A,A,MOD)
        k >>= 1

    return D

A=mat_power(a,K,MOD)
ans=0
for i in range(N):
    for j in range(N):
        ans=(ans+A[i][j])%MOD

print(ans)

Solution (Python):
import numpy as np

def matmul(A,B):
    N = len(A)
    C = np.eye(N, dtype=int)
    for n in range(N):
        C[n] = (A[n,:][:,None] * B % mod).sum(axis=0) % mod
    return C

def mat_power(A,k):
    N = len(A)
    if k == 0:
        return np.eye(N,dtype=int)
    X = mat_power(A,k//2)
    X = matmul(X,X)
    return matmul(A,X) if k&1 else X

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

print(mat_power(A,K).sum()%mod)

Solution (Python):
import numpy as np

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

def mat_mul(A, B):
    A1, A2 = A >> 15, A & (1 << 15) - 1
    B1, B2 = B >> 15, B & (1 << 15) - 1
    X = np.dot(A1, B1) % MOD
    Y = np.dot(A2, B2) % MOD
    Z = (np.dot(A1 + A2, B1 + B2) - X - Y) % MOD
    return ((X << 30) + (Z << 15) + Y) % MOD
 
def mat_power(A, N):
    P = np.eye(A.shape[0], dtype=np.int64)
    while N:
        if N & 1:
            P = mat_mul(P, A)
        A = mat_mul(A, A)
        N >>= 1
    return P

A2 = np.array(A)
print((np.sum(mat_power(A2,K))%MOD))

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

def matpow(vec,n):
    B=[[0]*len(vec) for i in range(len(vec))]
    for i in range(len(vec)):
        B[i][i]=1
    while(n>0):
        if n&1:
            B=mul(B,vec)
        vec=mul(vec,vec)
        n>>=1
    return B

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

A_k=matpow(A,K)

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

Solution (Python):






import sys
input = sys.stdin.readline

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]



def DD(arg): return defaultdict(arg)


def inv(n): return pow(n, MOD-2, MOD)


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count




import numpy as np
N,K = IL()
logK = int_log(K,2)+1
data = ILL(N)
A = np.matrix(data,dtype=np.object)

ruiA = [A]
for i in range(logK+1):
  ruiA.append((ruiA[-1]*ruiA[-1])%MOD)
  
now = np.eye(N,dtype=np.object)
for b in range(logK+2)[::-1]:
  if (K>>b)&1:
    now = ruiA[b]*now
    now %= MOD

print(int(np.sum(now)%MOD))

Solution (Python):






import sys
input = sys.stdin.readline

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]



def DD(arg): return defaultdict(arg)


def inv(n): return pow(n, MOD-2, MOD)


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count





N,K = IL()
data = ILL(N)

logK = int_log(K,2)+1
db = [[[0 for i in range(N)] for j in range(N)] for b in range(logK+1)] 
for i in range(N):
  db[0] = data

for b in range(logK):
  for j in range(N):
    for i in range(N):
      db[b+1][j][i] = sum(db[b][j][k]*db[b][k][i] for k in range(N))%MOD
import copy
now = [[0 for i in range(N)] for j in range(N)]
for i in range(N):
  now[i][i] = 1
for b in range(logK+1)[::-1]:
  if (K>>b)&1:
    prev = copy.deepcopy(now)
    for j in range(N):
      for i in range(N):
        now[j][i] = sum(prev[j][k]*db[b][k][i] for k in range(N))%MOD

print(sum(sum(l) for l in now)%MOD)

Solution (Python):
from sys import stdin
import sys
import math
from functools import reduce
import functools
import itertools
from collections import deque,Counter,defaultdict
from operator import mul
import copy


import heapq
sys.setrecursionlimit(10**6)
INF =  float("inf")
import bisect
import statistics
mod = 10**9+7

def matrix_mul(a,b):
    n = len(a)
    c = [[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):
                c[i][j] += (a[i][k]*b[k][j]) % mod

            c[i][j] = c[i][j] % mod

    return c

def mul_pow(x, k):
    n = len(x)
    c = [[0 for j in range(n)] for i in range(n)]

    for i in range(n):
        c[i][i] = 1

    while (k > 0):
        if k & 1: c = matrix_mul(c, x)
        x = matrix_mul(x, x)
        k = k >> 1
    return c

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

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

a = mul_pow(a, K)

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

print(ans % mod)

Solution (Python):
import sys
input = sys.stdin.readline
n, kk = map(int,input().split())
C = [list(map(int,input().split())) for i in range(n)]
mod = 10**9 + 7

l = 100
A = [[[0] * n for i in range(n)] for j in range(l)]

for i in range(n):
    for j in range(n):
        if C[i][j] == 1:
            A[0][i][j] = 1

for i in range(1, l):
    for j in range(n):
        for k in range(n):
            for m in range(n):
                A[i][j][m] = (A[i][j][m] + A[i-1][j][k] * A[i-1][k][m]) % mod


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

ind = 0
while kk > 0:
    if kk % 2 == 1:
        ANS_NEW = [[0] * n for i in range(n)]
        for j in range(n):
            for k in range(n):
                for m in range(n):
                    ANS_NEW[j][m] = (ANS_NEW[j][m] + ANS[j][k] * A[ind][k][m]) % mod

        for i in range(n):
            for j in range(n):
                ANS[i][j] = ANS_NEW[i][j]
    kk //= 2
    ind += 1

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

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor
from operator import mul
from functools import reduce

sys.setrecursionlimit(2147483647)
INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007

n, k = LI()
G = LIR(n)

def mat_mul(a, b):
    I, J, K = len(a), len(b[0]), len(b)
    c = [[0] * J for _ in range(I)]
    for i in range(I):
        for j in range(J):
            for k in range(K):
                c[i][j] += a[i][k] * b[k][j]
            c[i][j] %= mod
    return c

def mat_pow(x, n):
    y = [[0] * len(x) for _ in range(len(x))]
    for i in range(len(x)):
        y[i][i] = 1

    while n > 0:
        if n & 1:
            y = mat_mul(x, y)
        x = mat_mul(x, x)
        n >>= 1
    return y

ans = 0
GG = mat_pow(G, k)
for i in range(n):
    for j in range(n):
        ans += GG[i][j]

print(ans % mod)

Solution (Python):
import numpy as np

MOD = 10 ** 9 + 7


def pow_matrix_mod(x, n, mod=MOD):
    if not n:
        return np.eye(len(x), dtype=np.object)
    if n % 2 == 0:
        return pow_matrix_mod(x @ x % mod, n // 2) % mod
    else:
        return x @ pow_matrix_mod(x @ x % mod, (n - 1) // 2) % mod


N, K = map(int, input().split())
A = np.array([tuple(map(int, input().split())) for _ in range(N)], dtype=np.object)
print(pow_matrix_mod(A, K).sum() % MOD)

Solution (Python):
from sys import stdin, stdout, setrecursionlimit
from collections import deque, defaultdict, Counter
from heapq import heappush, heappop
from functools import lru_cache
import math

setrecursionlimit(10**6)
rl = lambda: stdin.readline()
rll = lambda: stdin.readline().split()
rli = lambda: map(int, stdin.readline().split())
rlf = lambda: map(float, stdin.readline().split())

INF, NINF = float('inf'), float('-inf')
MOD = 10**9 + 7

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

def mat_exp(M, k):
	n = len(M)
	ans = [[0 for __ in range(n)] for _ in range(n)]
	for i in range(n): ans[i][i] = 1
	while k:
		if k % 2 == 1: ans = mat_mul(ans, M)
		k //= 2
		M = mat_mul(M, M)
	return ans

def main():
	n, k = rli()
	A = []
	for _ in range(n):
		A.append(list(rli()))
	M = mat_exp(A, k)
	ans = 0
	for i in range(n):
		for j in range(n):
			ans += M[i][j]
			ans %= MOD
	print(ans)
	stdout.close()

if __name__ == "__main__":
	main()

Solution (Python):
def power(x,n,mod):
  digit = len(str(bin(n)))-2
  a = len(x)
  dp = [x]*digit
  for i in range(1,digit):
    dp[i] = dot(dp[i-1],dp[i-1],mod)
  ans = [[0]*a for _ in range(a)]
  for i in range(a):
    ans[i][i] = 1
  for i in range(digit):
    if 1<<i & n:
      ans = dot(dp[i],ans,mod)
  return ans

def dot(A,B,mod):
  p = len(A)
  q = len(A[0])
  r = len(B)
  s = len(B[0])
  if q!=r:
    print('掛け算できません')
    return
  ans = [[0]*s for _ in range(p)]
  for i in range(p):
    for j in range(s):
      for k in range(q):
        ans[i][j] += A[i][k]*B[k][j]
        ans[i][j] %= mod
  return ans

def solve():
  N, K = map(int, input().split())
  mat = [list(map(int, input().split())) for _ in range(N)]
  init = [[0]*N for _ in range(N)]
  for i in range(N):
    init[i][i] = 1
  mod = 10**9+7
  dp = dot(init,power(mat,K,mod),mod)
  ans = sum(map(sum,dp))%mod
  return ans
print(solve())

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

def calc(mat1,mat2):
  n=len(mat1)
  ret=[[0]*n for _ in range(n)]
  for i in range(n):
    for j in range(n):
      for k in range(n):
        ret[i][j]+=mat1[i][k]*mat2[k][j]
        ret[i][j]%=mod
  return ret

n,k=map(int,input().split())
tmp=[[0]*n for _ in range(n)]
for i in range(n):
  tmp[i][i]=1
matrix=[list(map(int,input().split())) for _ in range(n)]
for i in range(60):
  if k&(1<<i):
    tmp=calc(tmp,matrix)
  matrix=calc(matrix,matrix)  
ans=0
for i in range(n):
  for j in range(n):
    ans+=tmp[i][j]
    ans%=mod
print(ans)