/**
* created : 30.10.2024 12:36:13
* author : wzj33300
*/
#include <bits/stdc++.h>
using namespace std;
#ifdef DEBUG
#include <algo/debug.h>
#else
#define debug(...) 42
#define assert(...) 42
#endif
using ll = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
using db = long double; // or double, if TL is tight
using str = string; // yay python!
// pairs
using pi = pair<int, int>;
using pl = pair<ll, ll>;
using pd = pair<db, db>;
#define mp make_pair
#define fi first
#define se second
// ^ lol this makes everything look weird but I'll try it
template <class T>
using vc = vector<T>;
template <class T, size_t SZ>
using AR = array<T, SZ>;
using vi = vc<int>;
using vb = vc<bool>;
using vl = vc<ll>;
using vd = vc<db>;
using vs = vc<str>;
using vpi = vc<pi>;
using vpl = vc<pl>;
using vpd = vc<pd>;
// vectors
// oops size(x), rbegin(x), rend(x) need C++17
#define sz(x) int((x).size())
#define bg(x) begin(x)
#define all(x) bg(x), end(x)
#define rall(x) x.rbegin(), x.rend()
#define sor(x) sort(all(x))
#define rsz resize
#define ins insert
#define pb push_back
#define eb emplace_back
#define ft front()
#define bk back()
#define rep(i, n) for (int i = 0; i < (n); ++i)
#define rep1(i, n) for (int i = 1; i < (n); ++i)
#define rep1n(i, n) for (int i = 1; i <= (n); ++i)
#define repr(i, n) for (int i = (n) - 1; i >= 0; --i)
#define rep_subset(t, s) \
for (ll t = (s); t >= 0; t = (t == 0 ? -1 : (t - 1) & (s)))
#define lb lower_bound
#define ub upper_bound
template <class T>
int lwb(vc<T> &a, const T &b) {
return int(lb(all(a), b) - bg(a));
}
template <class T>
int upb(vc<T> &a, const T &b) {
return int(ub(all(a), b) - bg(a));
}
// const int MOD = 998244353; // 1e9+7;
const int Big = 1e9; // not too close to INT_MAX
const ll BIG = 1e18; // not too close to LLONG_MAX
const db PI = acos((db)-1);
const int dx[4]{1, 0, -1, 0}, dy[4]{0, 1, 0, -1}; // for every grid problem!!
mt19937 rng((uint32_t)chrono::steady_clock::now().time_since_epoch().count());
template <class T>
using pqg = priority_queue<T, vector<T>, greater<T>>;
int pct(int x) { return __builtin_popcount(x); }
int pct(u32 x) { return __builtin_popcount(x); }
int pct(ll x) { return __builtin_popcountll(x); }
int pct(u64 x) { return __builtin_popcountll(x); }
int popcnt_mod_2(int x) { return __builtin_parity(x); }
int popcnt_mod_2(u32 x) { return __builtin_parity(x); }
int popcnt_mod_2(ll x) { return __builtin_parityll(x); }
int popcnt_mod_2(u64 x) { return __builtin_parityll(x); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 1, 2)
int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(u32 x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
int topbit(u64 x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
// (0, 1, 2, 3, 4) -> (-1, 0, 1, 0, 2)
int lowbit(int x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(u32 x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(ll x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
int lowbit(u64 x) { return (x == 0 ? -1 : __builtin_ctzll(x)); }
template <typename T>
T floor(T a, T b) {
return a / b - (a % b && (a ^ b) < 0);
}
template <typename T>
T ceil(T x, T y) {
return floor(x + y - 1, y);
}
template <typename T>
T bmod(T x, T y) {
return x - y * floor(x, y);
}
template <typename T>
pair<T, T> divmod(T x, T y) {
T q = floor(x, y);
return {q, x - q * y};
}
template <class T>
bool ckmin(T &a, const T &b) {
return b < a ? a = b, 1 : 0;
} // set a = min(a,b)
template <class T>
bool ckmax(T &a, const T &b) {
return a < b ? a = b, 1 : 0;
} // set a = max(a,b)
template <class T, class U>
T fstTrue(T lo, T hi, U f) {
++hi;
assert(lo <= hi); // assuming f is increasing
while (lo < hi) { // find first index such that f is true
T mid = lo + (hi - lo) / 2;
f(mid) ? hi = mid : lo = mid + 1;
}
return lo;
}
template <class T, class U>
T lstTrue(T lo, T hi, U f) {
--lo;
assert(lo <= hi); // assuming f is decreasing
while (lo < hi) { // find first index such that f is true
T mid = lo + (hi - lo + 1) / 2;
f(mid) ? lo = mid : hi = mid - 1;
}
return lo;
}
int mod1, mod2;
int base;
int pw1[400 * 400 + 1], pw2[400 * 400 + 1];
struct Hash {
int x1 = 0, x2 = 0;
friend Hash &operator+=(Hash &me, int x) {
me.x1 += pw1[x];
me.x1 %= mod1;
me.x2 += pw2[x];
me.x2 %= mod2;
return me;
}
friend Hash operator+(const Hash &x, const Hash &y) {
return Hash((x.x1 + y.x1) % mod1, (x.x2 + y.x2) % mod2);
}
friend bool operator!=(const Hash &x, const Hash &y) {
return x.x1 != y.x2 || x.x2 != y.x2;
}
};
const int N = 400 * 400 * 20 * 2,
mod = 1e9 + 7;
int ls[N], rs[N], cnt, val[N], pw[400 * 400 + 1], c[N];
Hash sum[N];
void update(int &rt, int l, int r, int x) {
++cnt;
ls[cnt] = ls[rt], rs[cnt] = rs[rt], sum[cnt] = sum[rt], val[cnt] = val[rt], c[cnt] = c[rt];
rt = cnt;
c[rt]++;
val[rt] += pw[x];
val[rt] %= mod;
sum[rt] += x;
if (l == r) {
return;
}
int mid = l + r >> 1;
if (x <= mid)
update(ls[rt], l, mid, x);
else
update(rs[rt], mid + 1, r, x);
}
bool cmp(int x, int y, int l, int r) {
if (l == r) return c[x] > c[y];
int mid = l + r >> 1;
if (sum[rs[x]] != sum[rs[y]])
return cmp(rs[x], rs[y], mid + 1, r);
else
return cmp(ls[x], ls[y], l, mid);
}
using P = pair<int, int>;
P LS(P x) {
return P{ls[x.fi], ls[x.se]};
}
P RS(P x) {
return P{rs[x.fi], rs[x.se]};
}
Hash SUM(P x) {
return sum[x.fi] + sum[x.se];
}
int VAL(P x) {
return (val[x.fi] + val[x.se]) % mod;
}
int C(P x) {
return (c[x.fi] + c[x.se]) % mod;
}
bool cmp(P x, P y, int l, int r) {
if (l == r) return C(x) > C(y);
int mid = l + r >> 1;
if (SUM(RS(x)) != SUM(RS(y)))
return cmp(RS(x), RS(y), mid + 1, r);
else
return cmp(LS(x), LS(y), l, mid);
}
void _sol() {
int n, m;
cin >> n >> m;
pw[0] = 1;
for (int i = 1; i <= n * n; i++) {
pw[i] = 1ll * pw[i - 1] * n * n % mod;
}
vc<vi> a(n, vi(n));
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
cin >> a[i][j];
vc<vc<int>> f(n, vc<int>(n));
cnt = 0;
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++) {
if (i == 0 && j == 0) {
} else if (i == 0)
f[i][j] = f[i][j - 1];
else if (j == 0)
f[i][j] = f[i - 1][j];
else
f[i][j] = cmp(f[i - 1][j], f[i][j - 1], 1, n * n) ? f[i - 1][j] : f[i][j - 1];
update(f[i][j], 1, n * n, a[i][j]);
debug(val[f[i][j]]);
}
vc<vc<int>> g(n, vc<int>(n));
vc<vc<P>> pre(n, vc<P>(n)), suf(pre);
repr(i, n) repr(j, n) {
if (i == n - 1 && j == n - 1) {
} else if (i == n - 1)
g[i][j] = g[i][j + 1];
else if (j == n - 1)
g[i][j] = g[i + 1][j];
else
g[i][j] = cmp(g[i + 1][j], g[i][j + 1], 1, n * n) ? g[i + 1][j] : g[i][j + 1];
pre[i][j] = suf[i][j] = P{f[i][j], g[i][j]};
update(g[i][j], 1, n * n, a[i][j]);
}
debug(f, g);
for (int i = 0; i < n; i++) {
for (int j = 1; j < n; j++) {
pre[i][j] = cmp(pre[i][j], pre[i][j - 1], 1, n * n) ? pre[i][j] : pre[i][j - 1];
}
for (int j = n - 2; j >= 0; j--) {
suf[i][j] = cmp(suf[i][j], suf[i][j + 1], 1, n * n) ? suf[i][j] : suf[i][j + 1];
}
}
while (m--) {
int x1, x2, y1, y2;
cin >> x1 >> x2 >> y1 >> y2;
--x1, --x2, --y1, --y2;
P ans{0, 0};
// (x2, y1)
if (x2 < n - 1 && y1 > 0) {
ans = cmp(ans, pre[x2 + 1][y1 - 1], 1, n * n) ? ans : pre[x2 + 1][y1 - 1];
}
// (x1, y2)
if (x1 > 0 && y2 < n - 1) {
ans = cmp(ans, suf[x1 - 1][y2 + 1], 1, n * n) ? ans : suf[x1 - 1][y2 + 1];
}
cout << VAL(ans) << '\n';
}
}
// signed main() {
int main() {
// freopen(".in", "r",stdin);
// freopen(".out","w",stdout);
ios::sync_with_stdio(false);
cin.tie(0);
base = rng();
mod1 = 998244353;
mod2 = 1e9 + 7;
pw1[0] = pw2[0] = 1;
for (int i = 1; i <= 400 * 400; i++) {
pw1[i] = pw1[i - 1] * base % mod1;
pw2[i] = pw2[i - 1] * base % mod2;
}
int t;
cin >> t;
while (t--) {
_sol();
}
return 0;
}