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):
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={'
    
    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

read=sys.stdin.buffer.read;readline=sys.stdin.buffer.readline;input=lambda:sys.stdin.readline().rstrip()













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={'
    
    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

read=sys.stdin.buffer.read
readline=sys.stdin.buffer.readline
input=lambda: sys.stdin.readline().rstrip()















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)

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


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)

Solution (Python):

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()

Solution (Python):




import sys
import numpy as np

sys.setrecursionlimit(10**6)
input = sys.stdin.buffer.readline
INF = 10 ** 9 + 1  
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():
    
    N, K = map(int, input().split())
    X = np.int64(read().split())
    X = X.reshape((N, N))
    print(solve(N, K, X))



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



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":
        
        
        from my_module import solve  
    elif sys.argv[-1] == "-t":
        print("testing")
        _test()
        sys.exit()
    elif sys.argv[-1] != '-p' and len(sys.argv) == 2:
        
        input_as_file = open(sys.argv[1])
        input = input_as_file.buffer.readline
        read = input_as_file.buffer.read

    main()

Solution (Python):
class SquareMatrix():
    def __init__(self, n, mod=1000000007):
        self.n = n
        self.mat = [[0 for j in range(n)] for i in range(n)]
        self.mod = mod

    @staticmethod
    def id(n, mod=1000000007):
        res = SquareMatrix(n, mod)
        for i in range(n):
            res.mat[i][i] = 1
        return res

    @staticmethod
    def modinv(n, mod):
        assert n % mod != 0
        c0, c1 = n, mod
        a0, a1 = 1, 0
        b0, b1 = 0, 1
        while c1:
            a0, a1 = a1, a0 - c0 // c1 * a1
            b0, b1 = b1, b0 - c0 // c1 * b1
            c0, c1 = c1, c0 % c1
        return a0 % mod

    def set(self, arr):
        for i in range(self.n):
            for j in range(self.n):
                self.mat[i][j] = arr[i][j] % self.mod

    def operate(self, vec):
        assert len(vec) == self.n
        res = [0 for _ in range(self.n)]
        for i in range(self.n):
            for j in range(self.n):
                res[i] += self.mat[i][j] * vec[j]
                res[i] %= self.mod
        return res

    def add(self, other):
        assert other.n == self.n
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[i][j] + other.mat[i][j]
                res.mat[i][j] %= self.mod
        return res

    def subtract(self, other):
        assert other.n == self.n
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[i][j] - other.mat[i][j]
                res.mat[i][j] %= self.mod
        return res

    def times(self, k):
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[i][j] * k
                res.mat[i][j] %= self.mod
        return res

    def multiply(self, other):
        assert self.n == other.n
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                for k in range(self.n):
                    res.mat[i][j] += self.mat[i][k] * other.mat[k][j]
                    res.mat[i][j] %= self.mod
        return res

    def power(self, k):
        tmp = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                tmp.mat[i][j] = self.mat[i][j]
        res = SquareMatrix.id(self.n, self.mod)
        while k:
            if k & 1:
                res = res.multiply(tmp)
            tmp = tmp.multiply(tmp)
            k >>= 1
        return res

    def trace(self):
        res = 0
        for i in range(self.n):
            res += self.mat[i][i]
            res %= self.mod
        return res

    def determinant(self):
        res = 1
        tmp = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                tmp.mat[i][j] = self.mat[i][j]
        for j in range(self.n):
            if tmp.mat[j][j] == 0:
                for i in range(j + 1, self.n):
                    if tmp.mat[i][j] != 0:
                        idx = i
                        break
                else:
                    return 0
                for k in range(self.n):
                    tmp.mat[j][k], tmp.mat[idx][k] = tmp.mat[idx][k], tmp.mat[j][k]
                res *= -1
            inv = SquareMatrix.modinv(tmp.mat[j][j], self.mod)
            for i in range(j + 1, self.n):
                c = -inv * tmp.mat[i][j] % self.mod
                for k in range(self.n):
                    tmp.mat[i][k] += c * tmp.mat[j][k]
                    tmp.mat[i][k] %= self.mod
        for i in range(self.n):
            res *= tmp.mat[i][i]
            res %= self.mod
        return res

    def transpose(self):
        res = SquareMatrix(self.n, self.mod)
        for i in range(self.n):
            for j in range(self.n):
                res.mat[i][j] = self.mat[j][i]
        return res

    def inverse(self): 
        res = SquareMatrix.id(self.n, self.mod)
        tmp = SquareMatrix(self.n, self.mod)
        sgn = 1
        for i in range(self.n):
            for j in range(self.n):
                tmp.mat[i][j] = self.mat[i][j]
        for j in range(self.n):
            if tmp.mat[j][j] == 0:
                for i in range(j + 1, self.n):
                    if tmp.mat[i][j] != 0:
                        idx = i
                        break
                else:
                    return 0
                for k in range(self.n):
                    tmp.mat[j][k], tmp.mat[idx][k] = tmp.mat[idx][k], tmp.mat[j][k]
                    res.mat[j][k], res.mat[idx][k] = res.mat[idx][k], res.mat[j][k]
            inv = SquareMatrix.modinv(tmp.mat[j][j], self.mod)
            for k in range(self.n):
                tmp.mat[j][k] *= inv
                tmp.mat[j][k] %= self.mod
                res.mat[j][k] *= inv
                res.mat[j][k] %= self.mod
            for i in range(self.n):
                c = tmp.mat[i][j]
                for k in range(self.n):
                    if i == j:
                        continue
                    tmp.mat[i][k] -= tmp.mat[j][k] * c
                    tmp.mat[i][k] %= self.mod
                    res.mat[i][k] -= res.mat[j][k] * c
                    res.mat[i][k] %= self.mod
        return res

    def linear_equations(self, vec): 
        return self.inverse().operate(vec)

    def print(self):
        print(*self.mat, sep='\n')

MOD = 1000000007

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

M = SquareMatrix(N)
M.set(A)

P = M.power(K)

res = 0

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

print(res)

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



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])
    
    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


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)

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

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

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

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

Solution (Python):
from collections import Counter

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

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

    maxD = K.bit_length() - 1

    dp = [[Counter() for _ in range(N)] for _ in range(maxD+1)]
    for v, As in enumerate(Ass):
        v2s = []
        for v2, A in enumerate(As):
            if A == 1:
                v2s.append(v2)
        dp[0][v] = Counter(v2s)
    for d in range(1, maxD+1):
        for v in range(N):
            for v2, num2 in dp[d-1][v].items():
                for v3, num3 in dp[d-1][v2].items():
                    dp[d][v][v3] += num2 * num3
                    dp[d][v][v3] %= MOD

    anss = [1] * N
    for d in range(maxD+1):
        if K & (1<<d):
            ans2s = [0] * N
            for v in range(N):
                for v2, num in dp[d][v].items():
                    ans2s[v2] += anss[v] * num
                    ans2s[v2] %= MOD
            anss = ans2s

    print(sum(anss) % MOD)


solve()

Solution (Python):


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())

Solution (Python):


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())

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

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

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

print(ans.sum() % mod)