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

input = sys.stdin.readline


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


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

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


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

Solution (Python):
from numpy import*
M=10**9+7
m=lambda*z:map(int,raw_input().split())
n,k=m()
b=[1]*n
a=matrix(map(m,b),'O')
while k:b*=[1,a][k&1];a=a*a%M;k/=2
print sum(b)%M

Solution (Python):
from numpy import*
M=10**9+7
m=lambda*z:map(int,raw_input().split())
n,k=m()
b=[1]*n
a=matrix(map(m,b),'O')
while k:b=[b,b*a][k&1];a=a*a%M;k/=2
print sum(b)%M

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,raw_input().split())
n,k=m()
b=eye(n,n,0,'O')
a=matrix([m()for i in b],'O')
while k:b,a=[b,b*a][k&1],a*a%M;k/=2
print sum(b)%M

Solution (Python):
from numpy import*
M=10**9+7
m=lambda*z:map(int,raw_input().split())
n,k=m()
b=eye(n,n,0,'O')
a=matrix(map(m,b),'O')
while k:b=[b,b*a][k&1];a=a*a%M;k/=2
print sum(b)%M

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,raw_input().split())
n,k=m()
b=eye(n,n,0,'O')
a=matrix([m()for i in b],'O')
while k:b=[b,b*a][k&1];a=a*a%M;k/=2
print sum(b)%M

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,raw_input().split())
n,k=m()
b=eye(n,n,0,'O')
a=matrix([m()for i in b],'O')
while k:b=[b,b*a][k&1];a=a*a%M;k>>=1
print sum(b)%M

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
b=eye(n,n,0,'O')
a=matrix([list(m())for i in b],'O')
while k:b=[b,b*a][k&1];a=a*a%M;k>>=1
print(sum(b)%M)

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
a=matrix([list(m())for i in[0]*n],'O')
b=eye(n,n,0,'O')
while k:b=[b,b*a][k&1];a=a*a%M;k>>=1
print(sum(b)%M)

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
a=matrix([list(m())for i in[0]*n],'O')
b=eye(n,n,0,'O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(sum(b)%M)

Solution (Python):
from numpy import*
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
a=matrix([list(m())for i in[0]*n],'O')
b=eye(n,n,0,'O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy as P
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
a=P.matrix([list(m())for _ in[0]*n],dtype='O')
b=P.eye(n,n,0,'O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy as P
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
a=P.matrix([list(m())for i in[0]*n],'O')
b=P.eye(n,n,0,'O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy as P
M=10**9+7
n,k=map(int,input().split())
a=P.matrix([list(map(int,input().split()))for i in[0]*n],'O')
b=P.eye(n,n,0,'O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy as P
M=10**9+7
m=lambda:map(int,input().split())
n,k=m()
a=P.matrix([list(m())for _ in[0]*n],dtype='O')
b=P.eye(n,dtype='O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy as P
M=10**9+7
n,k=map(int,input().split())
a=P.matrix([list(map(int,input().split()))for i in[0]*n],dtype='O')
b=P.eye(n,dtype='O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy
M=10**9+7
n,k=map(int,input().split())
a=numpy.matrix([list(map(int,input().split()))for i in[0]*n],dtype='O')
b=numpy.eye(n,dtype='O')
while k:b=[b,b*a%M][k&1];a=a*a%M;k>>=1
print(b.sum()%M)

Solution (Python):
import numpy as P
M=10**9+7
n,k=map(int,input().split())
a=P.matrix([list(map(int,input().split()))for i in[0]*n],dtype='O')
b=P.eye(n,dtype='O')
while k:
 if k&1:b=b*a%M
 a=a*a%M
 k>>=1
print(b.sum()%M)

Solution (Python):
import numpy
M=10**9+7
n,k=map(int,raw_input().split())
a=numpy.matrix([map(int,raw_input().split())for i in[0]*n],dtype='O')
b=numpy.eye(n,dtype='O')
while k:
 if k&1:b=b*a%M
 a=a*a%M
 k/=2
print b.sum()%M

Solution (Python):
import numpy
M=10**9+7
n,k=map(int,input().split())
a=numpy.matrix([list(map(int,input().split()))for i in[0]*n],dtype='O')
b=numpy.eye(n,dtype='O')
while k:
 if k&1:b=b*a%M
 a=a*a%M
 k>>=1
print(b.sum()%M)

Solution (Python):

import sys
input = sys.stdin.readline

mod = int(1e9 + 7)

def f(N, K, A):
    m = [[0] * N for _ in range(N)]
    for i in range(N):
        m[i][i] = 1
    k = K
    a = A

    def mulmod(a, b):
        m = [[0] * N for _ in range(N)]
        for i, r, in enumerate(a):
            for j, c in enumerate(zip(*b)):
                m[i][j] = sum(x * y for x, y in zip(r, c)) % mod
        return(m)

    while k > 0:
        if k & 1:
            m = mulmod(m, a)
        a = mulmod(a, a)
        k >>= 1

    return(sum(sum(x) for x in m) % mod)

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

print(f(N, K, a))

Solution (Python):

import sys
input = sys.stdin.readline

mod = int(1e9 + 7)

def f(N, K, A):
    m = [[0] * N for _ in range(N)]
    for i in range(N):
        m[i][i] = 1
    k = K
    a = A

    def mulmod(a, b):
        m = [[0] * N for _ in range(N)]
        for i, r, in enumerate(a):
            for j, c in enumerate(zip(*b)):
                m[i][j] = sum(x * y for x, y in zip(r, c)) % mod
        return(m)

    while k > 0:
        if k & 1:
            m = mulmod(m, a)
        a = mulmod(a, a)
        k >>= 1

    return(sum(sum(x) for x in m) % mod)

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

print(f(N, K, a))

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

def MatMul(A,B):
    C=[[0]*len(B[0]) for _ in range(len(A))]
    for i in range(len(A)):
        for j in range(len(B[0])):
            for k in range(len(B)):
                C[i][j]+=A[i][k]*B[k][j]
                C[i][j]%=mod
    return C
def MatPow(A,n):
    l=len(A)
    B=[[0]*l for _ in range(l)]
    for i in range(l):
        B[i][i]=1
    while n:
        if n&1:
            B=MatMul(B,A)
        A=MatMul(A,A)
        n//=2
    return B

n,k=map(int,input().split())
A=[list(map(int,input().split())) for _ in range(n)]
P=MatPow(A,k)
ans=0
for p in P:
    ans+=sum(p)
    ans%=mod
print(ans)

Solution (Python):
n, k = map(int, input().split())
a = [list(map(int, input().split())) for _ in range(n)]

mod = 10**9 + 7

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

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

A = matpow(a, k)
ans = 0
for Ai in A:
    for Aij in Ai:
        ans += Aij
        ans %= mod

print(ans)

Solution (Python):


MOD = 1000000007


def mult(a, b):
    c = [[0] * len(a) for _ 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] += a[i][k] * b[k][j] % MOD
                c[i][j] %= MOD
    return c

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

    while n:
        if n % 2:
            p = mult(p, a)
        a = mult(a, a)
        n /= 2
    return p

n, k = map(int, raw_input().split())
graph = []
for i in range(n):
    graph.append(map(int, raw_input().split()))
graph = power(graph, k)
resp = 0
for r in graph:
    resp += sum(r)
    resp %= MOD
print(resp)

Solution (Python):

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

MOD = 10**9 + 7


dp = [ [0] * (N+1) for _ in range(N+1)]

S = [[1]*N]
KDIG = [0]*70

for i in range(1,61):
    if K % (1<<i) :
        KDIG[i-1] = 1
        K -= K%(1<<i)


def prod(A,B,k,l,m):
 
    C=[[None for i in range(m)] for j in range(k)]
 
    for i in range(k):
        for j in range(m):
            ANS=0
            for pl in range(l):
                ANS=(ANS+A[i][pl]*B[pl][j])%MOD
 
            C[i][j]=ANS
 
    return C



for i in range(60):
    if KDIG[i] == 1:
        S = prod(S,A,1,N,N)

    A = prod(A,A,N,N,N)
print(sum(S[0])%MOD)

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

S=[[1]*N]
KDIG=[0]*70

for i in range(1,61):
    if K%(1<<i):
        KDIG[i-1]=1
        K-=K%(1<<i)

def prod(A,B,k,l,m):

    C=[[None for i in range(m)] for j in range(k)]

    for i in range(k):
        for j in range(m):
            ANS=0
            for pl in range(l):
                ANS=(ANS+A[i][pl]*B[pl][j])%mod

            C[i][j]=ANS

    return C

for i in range(60):
    if KDIG[i]==1:
        S=prod(S,A,1,N,N)
    
    A=prod(A,A,N,N,N)

print(sum(S[0])%mod)

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

def mul(L,R):
    ret = [[0]*N for _ in range(N)]
    for i,row in enumerate(L):
        for j,col in enumerate(zip(*R)):
            ret[i][j] = sum(a*b for a,b in zip(row,col)) % MOD
    return ret

M = A
Z = [[int(i==j) for j in range(N)] for i in range(N)]
while K:
    if K%2:
        Z = mul(Z,M)
    M = mul(M,M)
    K //= 2

print(sum(sum(row) for row in Z) % MOD)

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


class Matrix():
    def __init__(self, A):
        self.A = A
        self.row = len(A)
        self.col = len(A[0])

    def __iter__(self):
        return self.A.__iter__()

    def __getitem__(self, i):
        return self.A.__getitem__(i)

    def __add__(self, B):
        aa = self.A
        bb = B.A
        return Matrix([[aa[i][j] + bb[i][j] for j in range(self.col)] for i in range(self.row)])

    def __sub__(self, B):
        aa = self.A
        bb = B.A
        return Matrix([[aa[i][j] - bb[i][j] for j in range(self.col)] for i in range(self.row)])

    def __mul__(self, B):
        bb = [[B.A[j][i] for j in range(B.row)] for i in range(B.col)]
        return Matrix([[sum([ak * bk for ak,bk in zip(ai,bj)]) % mod for bj in bb] for ai in self.A])

    def __truediv__(self, x):
        pass

    def pow(self, n):
        A = self
        r = Matrix([[0 if j != i else 1 for j in range(self.row)] for i in range(self.row)])
        while n > 0:
            if n % 2 == 1:
                r = r * A
            A = A * A
            n //= 2
        return r

    def lu(self):
        
        size = self.row
        T = copy.deepcopy(self.A)
        L = [[0]*size for _ in range(size)]
        U = [[0]*size for _ in range(size)]
        for i in range(size):
            for j in range(i,size):
                L[j][i] = T[j][i]
            for j in range(i,size):
                U[i][j] = T[i][j] / T[i][i]
            for j in range(i+1,size):
                for k in range(i+1,size):
                    T[j][k] -= L[j][i] * U[i][k]

        return Matrix(L),Matrix(U)

    def __str__(self):
        return self.A.__str__()


def main():
    n,k = LI()
    a = [LI() for _ in range(n)]
    A = Matrix(a)
    P = A.pow(k)

    return sum(map(sum,P.A)) % mod


print(main())

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


n, k = map(int, input().split())
dp = [list(map(int, input().split())) for _ in range(n)]
ans = [[0] * n for _ in range(n)]
for i in range(n):
    ans[i][i] = 1
MOD = 1000000007
while k:
    k, i = divmod(k, 2)
    if i:
        ans = dot(ans, dp)
    dp = dot(dp, dp)
print(sum(sum(a) % MOD for a in ans) % MOD)

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

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

R = [[i==j for i in range(N)] for j in range(N)]
while K:
    if K & 1:
        R = mul(R, A)
    A = mul(A, A)
    K >>= 1
print(sum(sum(r) for r in R) % MOD)

Solution (Python):
n,k = map(int,input().split())
g = [list(map(int,input().split())) for i in range(n)]
mod = 10**9+7
ans = 0
def mul(x,y):
  res = []
  for i in range(n):
    tmp = []
    for j in range(n):
      ko = 0
      for k in range(n):
        ko += x[i][k] * y[k][j]
        ko %=mod
      tmp.append(ko)
    res.append(tmp)
  return res
while k>0:
  if k &1 ==1:
    if type(ans) is int:
      ans = g
    else:
      ans = mul(g,ans)
  g = mul(g,g)
  k >>= 1
res = 0
for i in range(n):
  for j in range(n):
    res += ans[i][j]
    res %= mod
print(res)

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


def dot(X, Y):
    Z = [[0] * N for i in range(N)]
    for i in range(N):
        for j in range(N):
            t = 0
            for k in range(N):
                t += X[i][k] * Y[k][j]
                t %= mod
            Z[i][j] = t
    return Z


x = K
T = A
ans = [[0]*N for i in range(N)]
for i in range(N):
    ans[i][i]=1
while x > 0:
    if x % 2 == 1:
        ans = dot(ans, T)
    T = dot(T, T)
    
    x = x // 2
num = 0
for i in range(N):
    num += sum(ans[i])
    num %= mod
print(num)