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.



Sample input:
4 2
0 1 0 0
0 0 1 1
0 0 0 1
1 0 0 0


Sample output:
6

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

Solution code:
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()

Submission status: Accepted

Solution code:
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())

Submission status: Accepted

Solution code:
import sys,bisect,string,math,time,functools,random,fractions
from heapq import heappush,heappop,heapify
from collections import deque,defaultdict,Counter
from itertools import permutations,combinations,groupby
rep=range;R=range
def Golf():n,*t=map(int,open(0).read().split())
def I():return int(input())
def S_():return input()
def IS():return input().split()
def LS():return [i for i in input().split()]
def MI():return map(int,input().split())
def LI():return [int(i) for i in input().split()]
def LI_():return [int(i)-1 for i in input().split()]
def NI(n):return [int(input()) for i in range(n)]
def NI_(n):return [int(input())-1 for i in range(n)]
def StoLI():return [ord(i)-97 for i in input()]
def ItoS(n):return chr(n+97)
def LtoS(ls):return ''.join([chr(i+97) for i in ls])
def RA():return map(int,open(0).read().split())
def GI(V,E,ls=None,Directed=False,index=1):
    org_inp=[];g=[[] for i in range(V)]
    FromStdin=True if ls==None else False
    for i in range(E):
        if FromStdin:
            inp=LI()
            org_inp.append(inp)
        else:
            inp=ls[i]
        if len(inp)==2:
            a,b=inp;c=1
        else:
            a,b,c=inp
        if index==1:a-=1;b-=1
        aa=(a,c);bb=(b,c);g[a].append(bb)
        if not Directed:g[b].append(aa)
    return g,org_inp
def GGI(h,w,search=None,replacement_of_found='.',mp_def={'#':1,'.':0},boundary=1):
    #h,w,g,sg=GGI(h,w,search=['S','G'],replacement_of_found='.',mp_def={'#':1,'.':0},boundary=1) # sample usage
    mp=[boundary]*(w+2);found={}
    for i in R(h):
        s=input()
        for char in search:
            if char in s:
                found[char]=((i+1)*(w+2)+s.index(char)+1)
                mp_def[char]=mp_def[replacement_of_found]
        mp+=[boundary]+[mp_def[j] for j in s]+[boundary]
    mp+=[boundary]*(w+2)
    return h+2,w+2,mp,found
def TI(n):return GI(n,n-1)
def accum(ls):
    rt=[0]
    for i in ls:rt+=[rt[-1]+i]
    return rt
def bit_combination(n,base=2):
    rt=[]
    for tb in R(base**n):s=[tb//(base**bt)%base for bt in R(n)];rt+=[s]
    return rt
def gcd(x,y):
    if y==0:return x
    if x%y==0:return y
    while x%y!=0:x,y=y,x%y
    return y
def YN(x):print(['NO','YES'][x])
def Yn(x):print(['No','Yes'][x])
def show(*inp,end='\n'):
    if show_flg:print(*inp,end=end)

mo=10**9+7
inf=float('inf')
FourNb=[(-1,0),(1,0),(0,1),(0,-1)];EightNb=[(-1,0),(1,0),(0,1),(0,-1),(1,1),(-1,-1),(1,-1),(-1,1)];compas=dict(zip('WENS',FourNb));cursol=dict(zip('LRUD',FourNb))
l_alp=string.ascii_lowercase
#sys.setrecursionlimit(10**7)
read=sys.stdin.buffer.read;readline=sys.stdin.buffer.readline;input=lambda:sys.stdin.readline().rstrip()

########################################################################################################################################################################
# Verified by 
# https://atcoder.jp/contests/abc014/tasks/abc014_4
# https://atcoder.jp/contests/abc133/tasks/abc133_f

# import SparseTable
# initialize
# Tree() => n; (a,b);*(n-1)
# Tree(n) => (a,b);*(n-1)
# Tree(init=False); Tree.stdin(); => n; (a,b);*(n-1)
# Tree(init=False); Tree.listin(ls,index=0); => (a,b);*(n-1)

import sys,bisect,string,math,time,functools,random,fractions
from heapq import heappush,heappop,heapify
from collections import deque,defaultdict,Counter
from itertools import permutations,combinations,groupby
rep=range
def Golf():n,*t=map(int,open(0).read().split())
def I():return int(input())
def S_():return input()
def IS():return input().split()
def LS():return [i for i in input().split()]
def MI():return map(int,input().split())
def LI():return [int(i) for i in input().split()]
def LI_():return [int(i)-1 for i in input().split()]
def NI(n):return [int(input()) for i in range(n)]
def NI_(n):return [int(input())-1 for i in range(n)]
def StoLI():return [ord(i)-97 for i in input()]
def ItoS(n):return chr(n+97)
def LtoS(ls):return ''.join([chr(i+97) for i in ls])
def RA():return map(int,open(0).read().split())
def GI(V,E,ls=None,Directed=False,index=1):
    org_inp=[];g=[[] for i in range(V)]
    FromStdin=True if ls==None else False
    for i in range(E):
        if FromStdin:
            inp=LI()
            org_inp.append(inp)
        else:
            inp=ls[i]
        if len(inp)==2:
            a,b=inp;c=1
        else:
            a,b,c=inp
        if index==1:a-=1;b-=1
        aa=(a,c);bb=(b,c);g[a].append(bb)
        if not Directed:g[b].append(aa)
    return g,org_inp
def GGI(h,w,search=None,replacement_of_found='.',mp_def={'#':1,'.':0},boundary=1):
    #h,w,g,sg=GGI(h,w,search=['S','G'],replacement_of_found='.',mp_def={'#':1,'.':0},boundary=1) # sample usage
    mp=[boundary]*(w+2);found={}
    for i in range(h):
        s=input()
        for char in search:
            if char in s:
                found[char]=((i+1)*(w+2)+s.index(char)+1)
                mp_def[char]=mp_def[replacement_of_found]
        mp+=[boundary]+[mp_def[j] for j in s]+[boundary]
    mp+=[boundary]*(w+2)
    return h+2,w+2,mp,found
def TI(n):return GI(n,n-1)
def accum(ls):
    rt=[0]
    for i in ls:rt+=[rt[-1]+i]
    return rt
def bit_combination(n,base=2):
    rt=[]
    for tb in range(base**n):s=[tb//(base**bt)%base for bt in range(n)];rt+=[s]
    return rt
def gcd(x,y):
    if y==0:return x
    if x%y==0:return y
    while x%y!=0:x,y=y,x%y
    return y
def YN(x):print(['NO','YES'][x])
def Yn(x):print(['No','Yes'][x])
def show(*inp,end='\n'):
    if show_flg:print(*inp,end=end)

mo=10**9+7
inf=float('inf')
FourNb=[(-1,0),(1,0),(0,1),(0,-1)];EightNb=[(-1,0),(1,0),(0,1),(0,-1),(1,1),(-1,-1),(1,-1),(-1,1)];compas=dict(zip('WENS',FourNb));cursol=dict(zip('LRUD',FourNb))
l_alp=string.ascii_lowercase
#sys.setrecursionlimit(10**7)
read=sys.stdin.buffer.read
readline=sys.stdin.buffer.readline
input=lambda: sys.stdin.readline().rstrip()

## Verified by Yukicoder 1073
## https://yukicoder.me/problems/no/1073
##
## Matrix Class supporting operators +, -, *, %, +=, -=, *=, %=
## *, *= allows int/float/complex
## ** or pow(self,p,mod) for the size N*N matrix is implemented by Repeated squaring. O(N^3*log(p))
##
## Constructor: matrix(array), where array is 1D or 2D array. 1-dimensional array X is modified as 2D array of [X].
##
## methods
## T(): returns transposed matrix
## resize((n,m),fill=0): changes the matrix instance into the new shape (n * m), missing entries are filled with "fill" (default value is zero).


class matrix:
    class MulShapeError(Exception):
        "mult is not applicable between the two matrices given"
        pass

    def __init__(self,arr_input):
        if hasattr(arr_input[0],"__getitem__"):
            self.arr=arr_input
        else:
            self.arr=[arr_input]
        self.shape=(len(self.arr),len(self.arr[0]))

    def __getitem__(self,key):
        return self.arr[key]

    def __setitem__(self,key,value):
        self.arr[key]=value

    def __iter__(self):
        return iter(self.arr)
        
    def __add__(self,B):
        if type(B)!=matrix:
            return NotImplemented
        if B.shape!=self.shape:
            return NotImplemented
        rt=[[0]*self.shape[1] for i in range(self.shape[0])]
        for i in range(self.shape[0]):
            for j in range(self.shape[1]):
                rt[i][j]=self.arr[i][j]+B.arr[i][j]
        return matrix(rt)

    def __iadd__(self,B):
        return self.__add__(B)

    def __sub__(self,B):
        if type(B)!=matrix:
            return NotImplemented
        if B.shape!=self.shape:
            return NotImplemented
        rt=[[0]*self.shape[1] for i in range(self.shape[0])]
        for i in range(self.shape[0]):
            for j in range(self.shape[1]):
                rt[i][j]=self.arr[i][j]-B.arr[i][j]
        return matrix(rt)

    def __isub__(self,B):
        return self.__sub__(B)

    def __mul__(self,M):
        if type(M) in [int,float,complex]:
            M=matrix([[M*(i==j) for j in range(self.shape[1])] for i in range(self.shape[1])])
        if type(M)!=matrix:
            return NotImplemented
        if M.shape[0]!=self.shape[1]:
            raise matrix.MulShapeError("mult is not applicable between the matrices of shape "+str(self.shape)+" and "+str(M.shape))
        ra,ca=self.shape
        rb,cb=M.shape
        c=[[0]*cb for i in range(ra)]
        for i in range(ra):
            for j in range(cb):
                for k in range(ca):
                    c[i][j]+=self.arr[i][k]*M.arr[k][j]
        return matrix(c)

    def __imul__(self,M):
        return self.__mul__(M)

    def __rmul__(self,M):
        if type(M) in [int,float,complex]:
            M=matrix([[M*(i==j) for j in range(self.shape[1])] for i in range(self.shape[1])])
        if type(M)!=matrix:
            return NotImplemented
        if M.shape[0]!=self.shape[1]:
            raise matrix.MulShapeError("mult is not applicable between the matrix shape "+str(self.shape)+" and "+str(M.shape))
        ra,ca=M.shape
        rb,cb=self.shape
        c=[[0]*cb for i in range(ra)]
        for i in range(ra):
            for j in range(cb):
                for k in range(ca):
                    c[i][j]+=M.arr[i][k]*self.arr[k][j]
        return matrix(c)

    def __mod__(self,p):
        if type(p)!=int:
            return NotImplemented
        c=[[0]*self.shape[1] for i in range(self.shape[0])]
        for i in range(self.shape[0]):
            for j in range(self.shape[1]):
                c[i][j]=self.arr[i][j]%p
        return matrix(c)

    def __imod__(self,p):
        return self.__mod__(p)

    def __pow__(self,p,mod=10**9+7):
        if type(p)!=int or self.shape[0]!=self.shape[1]:
            return NotImplemented
        A=matrix(self.arr)
        R=matrix([[1*(i==j) for j in range(self.shape[0])] for i in range(self.shape[0])])
        while p>0:
            if p&1:
                R*=A
                R%=mod
            A*=A
            A%=mod
            p>>=1
        return R

    def __neg__(self):
        return self.__mul__(-1)

    def __str__(self):
        rt='['
        for i in self.arr:
            rt=rt+str(i)+",\n"
        return rt[:-2]+']'

    def T(self):
        rt=[[0]*self.shape[0] for i in range(self.shape[1])]
        for i in range(self.shape[0]):
            for j in range(self.shape[1]):
                rt[j][i]=self.arr[i][j]
        return matrix(rt)

    def resize(self,new_shape,fill=0):
        t_arr=[]
        for i in self.arr:
            t_arr+=i
        t_arr.reverse()
        n,m=new_shape
        self.shape=(n,m)
        self.arr=[[fill]*m for i in range(n)]

        for i in range(self.shape[0]):
            for j in range(self.shape[1]):
                if t_arr:
                    self.arr[i][j]=t_arr.pop()
        return

    def view(self):
        for i in self.arr:
            print(i)


show_flg=False
show_flg=True

ans=0

n,k=LI()
a=[]
for i in range(n):
    a+=LI(),
M=matrix(a)
M**=k

print(sum([sum(i)for i in M])%mo)

Submission status: Accepted

Solution code:
import sys

stdin = sys.stdin


def ni(): return int(ns())


def na(): return list(map(int, stdin.readline().split()))


def naa(N): return [na() for _ in range(N)]


def ns(): return stdin.readline().rstrip()  # ignore trailing spaces


def matrix_cal(A, B, mod):
    a = len(A)
    b = len(A[0])
    c = len(B)
    d = len(B[0])

    assert(b == c)
    ans = [[0] * d for _ in range(a)]
    for aa in range(a):
        for dd in range(d):
            v = 0
            for bb in range(b):
                v = (v + A[aa][bb] * B[bb][dd]) % mod
            ans[aa][dd] = v
    return ans


def matrix_pow(A, n, mod=0):
    if n == 1:
        return A
    a = matrix_pow(A, n//2, mod)
    ans = matrix_cal(a, a, mod)
    if n % 2 == 1:
        ans = matrix_cal(ans, A, mod)
    return ans


N, K = na()
A = naa(N)
mod = 10 ** 9 + 7

B = matrix_pow(A, K, mod)

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


print(ans)

Submission status: Accepted

Solution code:

def resolve():
    def mat_mul(a, b):
        H = len(a)
        W = len(b[0])
        K = len(a[0])

        C = [[0] * W for _ in range(H)]
        for i in range(H):
            for j in range(W):
                for k in range(K):
                    C[i][j] += a[i][k] * b[k][j]
                    C[i][j] %= MOD
        return C

    def mat_pow(a, p):
        n = len(a)
        ret = [[0] * n for _ in range(n)]
        for i in range(n):
            ret[i][i] = 1

        while p > 0:
            if p % 2 == 1:
                ret = mat_mul(ret, a)
            a = mat_mul(a, a)
            p >>= 1
        return ret

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

    ret = mat_pow(AB, K)
    ans = 0
    for i in range(N):
        for j in range(N):
            ans += ret[i][j]
            ans %= MOD

    print(ans)


if __name__ == "__main__":
    resolve()

Submission status: Accepted

Solution code:
#!/usr/bin/env python3

#from collections import defaultdict
#from heapq import heappush, heappop
import sys
import numpy as np

sys.setrecursionlimit(10**6)
input = sys.stdin.buffer.readline
INF = 10 ** 9 + 1  # sys.maxsize # float("inf")
MOD = 10 ** 9 + 7


def debug(*x):
    print(*x, file=sys.stderr)


def solve(N, K, X):
    def modmul(x, y):
        ret = np.zeros(x.shape, np.int64)
        for i in range(N):
            for j in range(N):
                v = x[i, :] * y[:, j]
                v %= MOD
                ret[i, j] = v.sum() % MOD
        return ret

    powK = np.eye(N, dtype=np.int64)
    while K:
        if K & 1:
            powK = modmul(powK, X)
        X = modmul(X, X)
        K //= 2
    return powK.sum() % MOD


def main():
    # parse input
    N, K = map(int, input().split())
    X = np.int64(read().split())
    X = X.reshape((N, N))
    print(solve(N, K, X))


# tests
T1 = """
4 2
0 1 0 0
0 0 1 1
0 0 0 1
1 0 0 0
"""


def test_T1():
    """
    >>> as_input(T1)
    >>> main()
    6
    """


T2 = """
3 3
0 1 0
1 0 1
0 0 0
"""


def test_T2():
    """
    >>> as_input(T2)
    >>> main()
    3
    """


T3 = """
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
"""


def test_T3():
    """
    >>> as_input(T3)
    >>> main()
    1
    """


T4 = """
1 1
0
"""


def test_T4():
    """
    >>> as_input(T4)
    >>> main()
    0
    """


T5 = """
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
"""


def test_T5():
    """
    >>> as_input(T5)
    >>> main()
    957538352
    """
# add tests above


def _test():
    import doctest
    doctest.testmod()


def as_input(s):
    "use in test, use given string as input file"
    import io
    global read, input
    f = io.StringIO(s.strip())

    def input():
        return bytes(f.readline(), "ascii")

    def read():
        return bytes(f.read(), "ascii")


USE_NUMBA = False
if (USE_NUMBA and sys.argv[-1] == 'ONLINE_JUDGE') or sys.argv[-1] == '-c':
    print("compiling")
    from numba.pycc import CC
    cc = CC('my_module')
    cc.export('solve', solve.__doc__.strip().split()[0])(solve)
    cc.compile()
    exit()
else:
    input = sys.stdin.buffer.readline
    read = sys.stdin.buffer.read

    if (USE_NUMBA and sys.argv[-1] != '-p') or sys.argv[-1] == "--numba":
        # -p: pure python mode
        # if not -p, import compiled module
        from my_module import solve  # pylint: disable=all
    elif sys.argv[-1] == "-t":
        print("testing")
        _test()
        sys.exit()
    elif sys.argv[-1] != '-p' and len(sys.argv) == 2:
        # input given as file
        input_as_file = open(sys.argv[1])
        input = input_as_file.buffer.readline
        read = input_as_file.buffer.read

    main()

Submission status: Accepted

Solution code:
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


n,k = map(int, input().split())
a = [None]*n
for i in range(n):
    a[i] = list(map(int, input().split()))
    
### 行列積 matrix multiplication
### 繰り返し2乗法

M = 10**9+7
def mul(a,b):
    """行列a,bの積
    aかbがNone のとき、単位行列として扱う
    """
    if a is None:
        return b
    elif b is None:
        return a
    n,m = len(a), len(a[0])
    k,l = len(b), len(b[0])
    # k==m
    out = [[0]*l for _ in range(n)]
    for i in range(n):
        for j in range(l):
            for p in range(m):
                out[i][j] += a[i][p]*b[p][j]
                out[i][j] %= M
    return out

# path^hを求める
i = 0
ans = None
tmp = a
while k>0:
    if k&1:
        ans = mul(ans, tmp)
    tmp = mul(tmp, tmp)
    k = k>>1
print(sum(sum(item)for item in ans) % M)

Submission status: Accepted

Solution code:
import numpy as np
MOD = 10**9+7

N, K = map(int, input().split())
A = np.zeros([N, N], dtype='int64')
for i in range(N):
    *A[i], = map(int, input().split())

def dot_MOD(x1, x2, MOD):
    ret = np.zeros((x1.shape[1],x2.shape[0]), dtype='int64')
    for i in range(x1.shape[0]):
        for j in range(x2.shape[1]):
            ret[i][j] = np.mod(x1[i,:] * (x2[:,j][np.newaxis:]), MOD, dtype='int64').sum()
    ret %= MOD
    return ret

ans = np.eye(N, dtype='int64')
while K:
    if K&1:
        ans = dot_MOD(ans,A,MOD)
    A = dot_MOD(A, A, MOD)
    K >>= 1
print(np.mod(ans.sum(), MOD, dtype='int64'))

Submission status: Accepted

Solution code:
import sys
import numpy as np
read = sys.stdin.buffer.read
input = sys.stdin.buffer.readline
MOD = 1_000_000_007
def dot_MOD(x1, x2, MOD):
    """行列の積とＭＯＤ

    ＭＯＤを取りつつ行列の積をもとめる。

    Args:
        x1 (2Dndarry of int): 被積の行列
        x2 (2Dndarry of int): 加積の行列

    Returns:
        ret (2Dndarry of int): 求まった行列

    Raises:
        例外の名前: 例外の説明 (例 : 引数が指定されていない場合に発生 )

    Yields:
        戻り値の型: 戻り値についての説明

    Examples:
        関数の使い方について記載

    Note:
        注意事項などを記載

    """
    ret = np.zeros((x1.shape[1],x2.shape[0]), dtype='int64')
    for i in range(x1.shape[0]):
        for j in range(x2.shape[1]):
            ret[i][j] = np.mod(x1[i,:] * (x2[:,j][np.newaxis:]), MOD, dtype='int64').sum()
    ret %= MOD
    return ret

def x2_MOD(x, MOD):
    """行列の自冪とＭＯＤ

    ＭＯＤを取りつつ行列の自冪をもとめる。

    Args:
        x1 (2Dndarry of int): 行列

    Returns:
        ret (2Dndarry of int): 求まった行列

    Raises:
        例外の名前: 例外の説明 (例 : 引数が指定されていない場合に発生 )

    Yields:
        戻り値の型: 戻り値についての説明

    Examples:
        関数の使い方について記載

    Note:
        注意事項などを記載

    """
    ret = np.empty(x.shape, dtype='int64')
    for i in range(x.shape[0]):
        for j in range(x.shape[0]):
            ret[i][j] = np.mod(x[i,:] * (x[:,j][np.newaxis:]), MOD, dtype='int64').sum()
    ret %= MOD
    return ret

N, K = map(int, input().split())
A = np.array(list(map(int, read().split())), dtype='int64').reshape((N, N))
ans = np.eye(N, dtype='int64')
while K:
    if K&1:
        ans = dot_MOD(ans,A,MOD)
    A = x2_MOD(A, MOD)
    K >>= 1
print(np.mod(ans.sum(), MOD, dtype='int64'))

Submission status: Accepted

Solution code:
from operator import mul

def solve():
    MOD = 10**9 + 7

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

    def getMatrixProduct(Ass, Bss, MOD):
        BssTr = map(list, zip(*Bss))
        Ax = len(Ass)
        By = len(Bss[0])
        ansss = [[0]*(By) for _ in range(Ax)]
        for y, Bs in enumerate(BssTr):
            for x, As in enumerate(Ass):
                ansss[x][y] = sum(map(mul, As, Bs)) % MOD
        return ansss
    def getMatrixPower(Ass, n, MOD):
        sizeA = len(Ass)
        ansss = [[0]*(sizeA) for _ in range(sizeA)]
        for i in range(sizeA):
            ansss[i][i] = 1
        while n:
            if n & 1:
                ansss = getMatrixProduct(ansss, Ass, MOD)
            Ass = getMatrixProduct(Ass, Ass, MOD)
            n //= 2
        return ansss

    Pss = getMatrixPower(Ass, K, MOD)

    ans = 0
    for Ps in Pss:
        for P in Ps:
            ans += P
            ans %= MOD

    print(ans)


solve()

Submission status: Accepted

Solution code:
import sys
readline = sys.stdin.buffer.readline
import numpy as np
sys.setrecursionlimit(10 ** 8)
INF = float('inf')
MOD = 10 ** 9 + 7


def main():
    N, K = map(int, readline().split())
    A = list(list(map(int, readline().split())) for _ in range(N))
    A = np.array(A)
    def dot(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 matrix_power(A, n):
        if n == 0:
            return np.eye(N, dtype=np.int64)
        X = matrix_power(A, n // 2)
        X = dot(X, X)
        return dot(A, X) if n & 1 else X

    ans = matrix_power(A,K)
    print(ans.sum() %MOD)


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:


def read_int():
    return int(input().strip())


def read_ints():
    return list(map(int, input().strip().split(' ')))


def merge(N, modulo, A, B):
    C = [[0]*N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            for j0 in range(N):
                C[i][j] = (C[i][j]+A[i][j0]*B[j0][j])%modulo
    return C


def solve():
    """
    OPT[k][i][j] - number of ways from i to j in K step
    OPT[1] = A

    OPT[k] = 
    """
    N, K = read_ints()
    A = []
    modulo = 10**9+7
    for _ in range(N):
        A.append(read_ints())
    def opt(k):
        if k == 1:
            return A
        half = opt(k//2)
        answer = merge(N, modulo, half, half)
        if k%2 == 1:
            answer = merge(N, modulo, answer, A)
        return answer
    return sum(c for row in opt(K) for c in row)%modulo


if __name__ == '__main__':
    print(solve())

Submission status: Accepted

Solution code:


def read_int():
    return int(input().strip())


def read_ints():
    return list(map(int, input().strip().split(' ')))


def solve():
    """
    OPT[k][i][j] - number of ways from i to j in K step
    OPT[1] = A

    OPT[k] = 
    """
    N, K = read_ints()
    A = []
    modulo = 10**9+7
    for _ in range(N):
        A.append(read_ints())
    def opt(k):
        if k == 1:
            return A
        half_A = opt(k//2)
        new_A = [[0]*N for _ in range(N)]
        for i in range(N):
            for j in range(N):
                for j0 in range(N):
                    new_A[i][j] = (new_A[i][j]+half_A[i][j0]*half_A[j0][j])%modulo
        if k%2 == 0:
            return new_A
        new_A2 = [[0]*N for _ in range(N)]
        for i in range(N):
            for j in range(N):
                for j0 in range(N):
                    new_A2[i][j] = (new_A2[i][j]+new_A[i][j0]*A[j0][j])%modulo
        return new_A2
    return sum(c for row in opt(K) for c in row)%modulo



if __name__ == '__main__':
    print(solve())

Submission status: Accepted

Solution code:
import numpy as np
MOD = 10 ** 9 + 7


def mult(N, 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, k):
    if k == 0:
        return np.eye(N, dtype=np.int64)
    X = mat_power(A, N, k//2)
    X = mult(N, X, X)
    return mult(N, 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)], np.int64)
print(mat_power(a, n, k).sum() % MOD)

Submission status: Accepted

Solution code:
def pow(mat, k, mod):
    n = len(mat)
    if k == 0:
        return [[1 if i == j else 0 for i in range(n)] for j in range(n)]
    if k % 2 == 0:
        sqrt_mat = pow(mat, k // 2, mod)
        # calculate 'sqrt_mat * sqrt_mat'
        ret = [[0 for _ in range(n)] for _ in range(n)]
        for i in range(n):
            for j in range(n):
                v = 0
                for k in range(n):
                    v += sqrt_mat[i][k] * sqrt_mat[k][j]
                    v %= mod
                ret[i][j] = v
        return ret
    else:
        powmat = pow(mat, k - 1, mod)
        # calculate 'mat * powmat'
        ret = [[0 for _ in range(n)] for _ in range(n)]
        for i in range(n):
            for j in range(n):
                v = 0
                for k in range(n):
                    v += mat[i][k] * powmat[k][j]
                    v %= mod
                ret[i][j] = v
        return ret


def main():
    MOD = 10 ** 9 + 7
    N, K = list(map(int, input().split(' ')))
    A = [list(map(int, input().split(' '))) for _ in range(N)]
    D = pow(A, K, MOD)
    ans = 0
    for i in range(N):
        for j in range(N):
            ans += D[i][j]
            ans %= MOD
    print(ans)


if __name__ == '__main__':
    main()

Submission status: Accepted

Solution code:
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:
    K, b = divmod(K, 2)
    if b:
        ans = ans.dot(A) % mod
    A = A.dot(A) % mod

print(ans.sum() % mod)

Submission status: Accepted

Solution code:
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)