思路:

我们可以预处理出一个矩阵\(P\),其中\(P[i][j]\)表示从\(i\)到\(j\)的道路数目.

这样我们大概可以搞出从\(u\)到\(v\)正好是\(d\)条路径的方案数.如果是不大于\(d\)条路径,我们就可以加一个计数器,每次累加能够到\(v\)的方案数.

这样就是一个\((N+1)*(N+1)\)的转移矩阵.如果我们每次用这个转移矩阵乘的话就要爆掉了.

发现矩阵\(P\)的算法非常奇葩.我们可以将其分解成一个\((N+1)*(K+1)\)的矩阵\(out\)与一个\((K+1)*(N+1)\)的矩阵\(in\)的乘积.

于是\(P^d=(out*in)^d\).

我们发现\(P^d=(out*in)^d=out*(in*out)^{d-1}*in\).注意到\(in*out\)是一个\((K+1)*(K+1)\)的矩阵,算起来毫无压力.

我们的答案矩阵是一个\(1*(N+1)\)的矩阵,容易发现这些东西乘到一起也只是\(O(NK^2)\)级别的.

因此我们的复杂度就是\(O(Q(K^3logd+NK^2))\).

(注意本题卡常数)

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
using namespace std;
static const int mod=(1e9)+7;
inline void inc(int&x,int y){if((x+=y)>=mod)x-=mod;}
struct Matrix{
    int d[1010][1010],w,h;
    Matrix(){}
    Matrix(int _w,int _h):w(_w),h(_h){memset(d,0,sizeof d);}
    void output(){
        printf("w=%d h=%d\n",w,h);
        for(int i=1;i<=w;++i){
            for(int j=1;j<=h;++j)printf("%d ",d[i][j]);puts("");
        }
        puts("");
    }
};
long long calc=0;
Matrix t,res,ret;
Matrix mul(const Matrix&A,const Matrix&B){
    res.w=A.w,res.h=B.h;
    for(int i=1;i<=res.w;++i)for(int j=1;j<=res.h;++j){
        res.d[i][j]=0;for(int k=1;k<=A.h;++k)++calc,inc(res.d[i][j],(long long)A.d[i][k]*B.d[k][j]%mod);
    }
    return res;
}
Matrix Pow(const Matrix&x,int y){//for w=h
    t=ret=x;y--;for(;y;y>>=1,t=mul(t,t))if(y&1)ret=mul(ret,t);return ret;
}
#define N 1010
#define M 21
int n,m,in[N][M],out[N][M];
Matrix Solve1,Solve2,Solve3,ans;
int main(){
    #ifndef ONLINE_JUDGE
    freopen("tt.in","r",stdin);
    #endif
    int n,m;scanf("%d%d",&n,&m);
    register int i,j;
    for(i=1;i<=n;++i){for(j=1;j<=m;++j)scanf("%d",&out[i][j]);for(j=1;j<=m;++j)scanf("%d",&in[i][j]);}
    int q,u,v,d;
    scanf("%d",&q);
    while(q--){
        calc=0;//debug
        scanf("%d%d%d",&u,&v,&d);
        if(d==0)printf("%d\n",u==v?1:0);
        else if(d==1){
            int ret=0;for(int i=1;i<=m;++i)ret=(ret+(long long)out[u][i]*in[v][i]%mod)%mod;
            printf("%d\n",(ret+(u==v?1:0))%mod);
        }
        else{
            Solve1=Matrix(n+1,m+1),Solve2=Matrix(m+1,n+1);
            for(i=1;i<=n;++i)for(j=1;j<=m;++j)Solve1.d[i][j]=out[i][j];
            Solve1.d[n+1][m+1]=1;
            for(i=1;i<=n;++i)for(j=1;j<=m;++j)Solve2.d[j][i]=in[i][j];
            for(j=1;j<=m;++j)Solve2.d[j][n+1]=in[v][j];
            Solve2.d[m+1][n+1]=1;
            //Solve1.output();
            //Solve2.output();
            Solve3=mul(Solve2,Solve1);//Solve3.output();
            Solve3=Pow(Solve3,d-1);
            ans=Matrix(1,n+1);ans.d[1][u]=1;
            ans=mul(ans,Solve1),ans=mul(ans,Solve3),ans=mul(ans,Solve2);
            printf("%d\n",(ans.d[1][n+1]+(u==v?1:0))%mod);
        }
        //printf("%I64d\n",calc);
    }
    return 0;
}

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
using namespace std;
static const int mod=(1e9)+7;
inline void inc(int&x,int y){if((x+=y)>=mod)x-=mod;}
#define N 1010
#define M 22
int n,m,in[N][M],out[N][M];
int Solve1[N][M],Solve2[M][N],Solve3[M][M],ans[N][N],reg[N][N],t[M][M],one[M][M];
inline void work(int d){//solve Solve3^d
    register int i,j,k;
    memset(one,0,sizeof one);for(i=1;i<=m+1;++i)one[i][i]=1;
    for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j)t[i][j]=Solve3[i][j];
    while(d){
        if(d&1){
            for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j){
                reg[i][j]=0;for(k=1;k<=m+1;++k)inc(reg[i][j],(long long)one[i][k]*t[k][j]%mod);
            }
            for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j)one[i][j]=reg[i][j];
        }
        d>>=1;
        for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j){
            reg[i][j]=0;for(k=1;k<=m+1;++k)inc(reg[i][j],(long long)t[i][k]*t[k][j]%mod);
        }
        for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j)t[i][j]=reg[i][j];
    }
    for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j)Solve3[i][j]=one[i][j];
}
void pr1(){
    puts("Solve1");
    for(int i=1;i<=n+1;++i){
        for(int j=1;j<=m+1;++j)printf("%9d ",Solve1[i][j]);puts("");
    }
}
void pr2(){
    puts("Solve2");
    for(int i=1;i<=m+1;++i){
        for(int j=1;j<=n+1;++j)printf("%9d ",Solve2[i][j]);puts("");
    }
}
void pr3(){
    puts("Solve3");
    for(int i=1;i<=m+1;++i){
        for(int j=1;j<=m+1;++j)printf("%9d ",Solve3[i][j]);puts("");
    }
}
int main(){
    #ifndef ONLINE_JUDGE
    freopen("tt.in","r",stdin);
    #endif
    scanf("%d%d",&n,&m);
    register int i,j,k;
    for(i=1;i<=n;++i){for(j=1;j<=m;++j)scanf("%d",&out[i][j]);for(j=1;j<=m;++j)scanf("%d",&in[i][j]);}
    int q,u,v,d;
    scanf("%d",&q);
    while(q--){
        scanf("%d%d%d",&u,&v,&d);
        if(d==0)printf("%d\n",u==v?1:0);
        else if(d==1){
            int ret=0;for(int i=1;i<=m;++i)ret=(ret+(long long)out[u][i]*in[v][i]%mod)%mod;
            printf("%d\n",(ret+(u==v?1:0))%mod);
        }
        else{
            memset(Solve1,0,sizeof Solve1);
            for(i=1;i<=n;++i)for(j=1;j<=m;++j)Solve1[i][j]=out[i][j];
            Solve1[n+1][m+1]=1;
            memset(Solve2,0,sizeof Solve2);
            for(i=1;i<=n;++i)for(j=1;j<=m;++j)Solve2[j][i]=in[i][j];
            for(j=1;j<=m;++j)Solve2[j][n+1]=in[v][j];
            Solve2[m+1][n+1]=1;
            
            for(i=1;i<=m+1;++i)for(j=1;j<=m+1;++j){
                Solve3[i][j]=0;for(k=1;k<=n+1;++k)inc(Solve3[i][j],(long long)Solve2[i][k]*Solve1[k][j]%mod);
            }
            
            work(d-1);
            
            for(j=1;j<=n+1;++j)ans[1][j]=0;ans[1][u]=1;
            
            for(i=1;i<=1;++i)for(j=1;j<=m+1;++j){
                reg[i][j]=0;for(k=1;k<=n+1;++k)inc(reg[i][j],(long long)ans[i][k]*Solve1[k][j]%mod);
            }for(i=1;i<=1;++i)for(j=1;j<=m+1;++j)ans[i][j]=reg[i][j];
            
            for(i=1;i<=1;++i)for(j=1;j<=m+1;++j){
                reg[i][j]=0;for(k=1;k<=m+1;++k)inc(reg[i][j],(long long)ans[i][k]*Solve3[k][j]%mod);
            }for(i=1;i<=1;++i)for(j=1;j<=m+1;++j)ans[i][j]=reg[i][j];
            
            for(i=1;i<=1;++i)for(j=1;j<=n+1;++j){
                reg[i][j]=0;for(k=1;k<=m+1;++k)inc(reg[i][j],(long long)ans[i][k]*Solve2[k][j]%mod);
            }for(i=1;i<=1;++i)for(j=1;j<=n+1;++j)ans[i][j]=reg[i][j];
            
            printf("%d\n",(ans[1][n+1]+(u==v?1:0))%mod);
        }
    }
    return 0;
}