[ again, this is for my personal usage. USACO solution is a little bit too difficult to understand, at least for me. The explanation below is easier to understand. ]

problem statement: http://usaco.org/index.php?page=viewproblem2&cpid=1069

dp(a,b,c) : number of ways from a -> b, the highest button used is upto c

• first, after we press the button, we can stay at b, so dp(a,b,c) += dp(a,b,c-1)

• or we go a -> r -> b with max button of c. In this case, between a->r, we find a point x with max button c-1. Then from x -> r the worse you can do is press button c between r -> b, since c-1 is reset and available, the worse we can do is press c-1 to a point y, then y->b use max button of c-1 then define dp1(a,r,c-1) is sum over all x, where x->r, dp(a,x,c-1) dp2(r,b,c-1) is sum over all y, where y->r, dp(y,b,c-1)

• then dp(a,b,c) += sum over r, dp1 * dp2

Then we put constrain on (bs,s) and (bt,t) use dp[2][2](a,b,c), 0/1 is whether s/t constrain is used or not Then we put constrain on (bs,s) and (bt,t) use dp[2][2](a,b,c), 0/1 is whether s/t constrain is used or not

• doing the same, a -> b is sum over r of a -> r -> b

• but a->r is pinned at the start, and r->b is pinned at the end

/*
ID: USACO_template
LANG: C++
PROG: USACO
*/
#include <iostream>  //cin , cout
#include <fstream>   //fin, fout
#include <stdio.h>   // scanf , pringf
#include <cstdio>
#include <algorithm> // sort , stuff
#include <stack>     // stacks
#include <queue>     // queues
#include <map>
#include <string.h>
#include <set>

using namespace std;

typedef pair<int, int>          pii;
typedef vector<int>             vi;     /// adjlist without weight
typedef vector<pii>             vii;    /// adjlist with weight
typedef vector<pair<int,pii>>   vpip;   /// edge with weight
typedef long long               ll;

#define mp  make_pair
#define fst first
#define snd second
#define pb  push_back
#define sz(x) (int)(x).size()
#define trav(u, adj_v) for (auto& u: adj_v)

const int MOD = 1e9+7;  // 998244353;
const int MX  = 2e5+5;   //
const ll  INF = 1e18;    //

#define MAXV 63
#define MAXE 100007

bool debug;

int N, K, Q;
ll adj[MAXV][MAXV];

/**
dp(a,b,c) : number of ways from a -> b, the highest button used is upto c
- first, after we press the button, we can stay at b, so dp(a,b,c) += dp(a,b,c-1)
- or we go a -> r -> b with max button of c. In this case,
    between a->r, we find a point x with max button c-1. Then from x -> r the worse you can do is press button c
    between r -> b, since c-1 is reset and available, the worse we can do is press c-1 to a point y, then y->b use max button of c-1
  then define
    dp1(a,r,c-1) is sum over all x, where x->r, dp(a,x,c-1)
    dp2(r,b,c-1) is sum over all y, where y->r, dp(y,b,c-1)
- then dp(a,b,c) += sum over r, dp1 * dp2

Then we put constrain on (bs,s) and (bt,t)
use dp[2][2](a,b,c), 0/1 is whether s/t constrain is used or not

*/
ll dp[2][2][MAXV][MAXV][MAXV];
ll dp1[2][MAXV][MAXV][MAXV], dp2[2][MAXV][MAXV][MAXV];

void InitDP() {
    for(int c = 1; c <=K; c++) {
        // calculate dp1 and dp2
        for(int r = 1; r<=N; r++) {
            for(int a = 1; a <=N; a++) {
                for(int x = 1; x<=N; x++) {
                    if(adj[x][r])
                        dp1[0][a][r][c] = (dp1[0][a][r][c] + dp[0][0][a][x][c-1]*1LL) % MOD;
                }
            }
            dp1[0][r][r][c] = ( dp1[0][r][r][c] + 1LL) % MOD;

            for(int y = 1; y <=N; y++) {
                for(int b=1; b<=N; b++) {
                    if(adj[r][y])
                        dp2[0][r][b][c] = (dp2[0][r][b][c] + 1LL*dp[0][0][y][b][c-1]) % MOD;
                }
            }
            dp2[0][r][r][c] = ( dp2[0][r][r][c] + 1LL) % MOD;
        }
        // calculate dp
        for(int a=1; a<=N; a++) {
            for(int b=1; b<=N; b++) {
                dp[0][0][a][b][c] = dp[0][0][a][b][c-1] % MOD;
                for(int r=1; r <=N; r++) {
                    dp[0][0][a][b][c] += ( dp1[0][a][r][c] * dp2[0][r][b][c] ) % MOD;
                    dp[0][0][a][b][c] %= MOD;
                }
            }
        }
    }

}

/**
Then we put constrain on (bs,s) and (bt,t)
use dp[2][2](a,b,c), 0/1 is whether s/t constrain is used or not
- doing the same, a -> b is sum over r of a -> r -> b
- but a->r is pinned at the start, and r->b is pinned at the end

*/
void initZero() {
    for(int i=0; i<2; i++) for(int j=0; j<2; j++) {
        if(i==0 && j==0) continue;
        memset(dp[i][j], 0LL, sizeof(dp[i][j]));
    }
    memset(dp1[1], 0LL, sizeof(dp1[1]));
    memset(dp2[1], 0LL, sizeof(dp2[1]));
}
ll calcDP(int bs, int s, int bt, int t) {
    initZero();
    for(int c=1; c<=K; c++) {
        for(int r=1; r<=N; r++) {
            for(int x=1; x<=N; x++) {
                if(adj[x][r])
                    dp1[1][s][r][c] = (dp1[1][s][r][c] + dp[1][0][s][x][c-1]*1LL) % MOD;
            }
            for(int y=1; y<=N; y++) {
                if(adj[r][y])
                    dp2[1][r][t][c] = (dp2[1][r][t][c] + 1LL*dp[0][1][y][t][c-1]) % MOD;
            }
        }

        if(c==bs) { dp1[1][s][s][c]++; dp1[1][s][s][c] %= MOD; }
        if(c==bt) { dp2[1][t][t][c]++; dp2[1][t][t][c] %= MOD; }

        dp[1][1][s][t][c] = (dp[1][1][s][t][c] + dp[1][1][s][t][c-1]) % MOD;
        for(int r=1; r<=N; r++) {
            dp[1][1][s][t][c] += ( dp1[1][s][r][c] * dp2[1][r][t][c] ) %MOD; dp[1][1][s][t][c] %= MOD;
        }

        // update for next iteration
        for(int r=1; r<=N; r++) {
            dp[1][0][s][r][c] += dp[1][0][s][r][c-1]; dp[1][0][s][r][c] %= MOD;
            dp[0][1][r][t][c] += dp[0][1][r][t][c-1]; dp[0][1][r][t][c] %= MOD;
            for(int r1=1; r1<=N; r1++) {
                dp[1][0][s][r][c] += (dp1[1][s][r1][c] * dp2[0][r1][r][c])%MOD; dp[1][0][s][r][c] %= MOD;
                dp[0][1][r][t][c] += (dp1[0][r][r1][c] * dp2[1][r1][t][c])%MOD; dp[0][1][r][t][c] %= MOD;
            }
        }

    }

	return dp[1][1][s][t][K];
}

int main() {
    debug = true;
    ios_base::sync_with_stdio(false); cin.tie(0);

    int i, j;

    cin >> N >> K >> Q;
    for(i=1; i<=N; i++)
        for(j=1;j<=N; j++) { char c; cin >> c; adj[i][j] = c - '0'; }

    InitDP();

    for(i=1; i<=Q; i++) {
        int bs, s, bt, t; cin >> bs >> s >> bt >> t;
        cout << calcDP(bs, s, bt, t) << endl;
    }

}