[ again, for a personal record. USACO explanation is 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;
}
}