题目概括：

P3474 [POI 2008] KUP-Plot purchase

给出两个正整数$k$与$n$， 以及一个$n\cdot n$的矩阵，请求出一个子矩阵使得其各项之和$\in[k,2k]$

如果您还不熟练求极大（最大）子矩阵的做法，建议您先去看如下几道题：

P4147 玉蟾宫

P1578 [WC2002] 奶牛浴场

以及其前置单调站相关知识：

CF817D Imbalanced Array

P2422 良好的感觉

P5788 【模板】单调栈

思路分析：

（1）.我们考虑先对所有数字扫一遍，对于任意一个$a_{i,j}$，总有如下三种情况

• $a_{i,j}<k$

• $k\le a_{i,j}\le 2k$ 即 $a_{i,j}\in[k,2k]$

• $a_{i,j}>2k$

对于情况1我们不做处理，对于情况2，我们的$a_{i,j}$就是一组解，这里需注意输出先列后行，因此输出应为

j i j i

对于情况3，我们将它们标记为“坏点”，因为如果我们选择的矩形纳入了它们中的任何一个，矩形的和一定大于 $2k$ 不满足题目条件，因此我们在之后的操作中一定不能选取它们

（2）.在极大子矩阵的和大于等于 $k$ 时，其中一定含有答案矩阵（极大子矩阵指的是不包含“坏点”的矩阵中矩阵之和最大的那个矩阵）

证明：

若极大子矩阵的和 $S$ 满足 $S \in [2,2k]$ ，其本身就是答案

否则，我们记极大子矩阵中的第 $i$ 行为 $M_i$

那么对于每个 $M_i$ ，有如下三种情况

• $M_i < k$

• $k \le M_i \le 2k$ 即 $M_i \in [k,2k]$

• $M_i > 2k$

对于情况2，那么它本身就是答案，输出即可，对于情况3，你还记得吗，在上文中我们提到了 $a_{i,j}$ 的处理，显然在处理后 $a_{i,j} < k$ ，于是我们可以遍历 $M_i$ ，来求得所需的答案，具体而言：

int sum = 0;
for (int j = s.second; j <= e.second; j++) {
	if (sum + a[i][j] > 2 * k) {
		cout << s.second << " " << i << " " << j - 1 << " " << i;
		return 0;
	}
	sum += a[i][j];
}

代码中的第3行是说若 $sum + a_{i,j}$ 即 $M_i$ 的前缀和大于 $2k$ 那么在没有加 $a_{i,j}$ 之前前缀和一定 $\in [k,2k]$ ，所以我们输出答案（其中 $s$ 表示极大矩阵的左上角，$e$ 表示极大矩阵的右下角）

记前 $j$ 项的前缀和为 $sum_j$ ，证明：

$$sum_j \ge 2k \\ M_{i,j} < k \\ sum_j = sum_{j-1}+M_{i,j} \\ sum_{j-1}=sum_j-M_{i,j} \ge sum_j-k \ge 2k-k \ge k \\ sum_{j-1} \le sum_j = 2k \\ k \le sum_{j-1} \le 2k$$

证毕

而对于情况1，类似于我们上述的处理，只要我们让极大子矩阵的和减掉 $M_i$ 总能够得到一个合法的矩阵，证明类比情况2

（3）.我们可以用直方图+单调站或者悬线法来求出极大子矩阵（本文使用直方图+单调站，其实这两个核心思想是一样的，实现略有不同）然后利用（2）中所说的方式求得答案，如果极大子矩阵的和小于 $k$ 输出

NIE

具体实现流程：

1.一些定义（其实开到2000就够，注意cnt一定要开long long）

int n, a[2001][2001], lmin[2001][2001], rmin[2001][2001], h[2001][2001];
long long cnt[2001][2001];

2.初始化，先扫一遍

int k;
cin >> k >> n;
for (int i = 1; i <= n; i++)
	for (int j = 1; j <= n; j++) {
		cin >> a[i][j];
		if (a[i][j] >= k && a[i][j] <= 2 * k) {
			cout << j << " " << i << " " << j << " " << i;
			return 0;
		}
		cnt[i][j] = cnt[i - 1][j] + cnt[i][j - 1] - cnt[i - 1][j - 1] + a[i][j];
	}

3.建立直方图

for (int i = 1; i <= n; i++)
	for (int j = 1; j <= n; j++)
		if (a[i][j] <= 2 * k)
			h[i][j] = h[i - 1][j] + 1;

3.求解每行的 $l$ (左侧第一个小于它的)与 $r$(右侧第一个小于它的)

void solve(int line) {
	stack<int> l, r;
	for (int i = 1; i <= n; i++) {
		while (!l.empty() && h[line][l.top()] >= h[line][i])
			l.pop();
		if (l.empty())
			lmin[line][i] = 0;
		else
			lmin[line][i] = l.top();
		l.push(i);
	}
	for (int i = n; i >= 1; i--) {
		while (!r.empty() && h[line][r.top()] >= h[line][i])
			r.pop();
		if (r.empty())
			rmin[line][i] = n + 1;
		else
			rmin[line][i] = r.top();
		r.push(i);
	}
}
...
for (int i = 1; i <= n; i++)
	solve(i);

4.寻找极大子矩阵

pair<int, int> s, e;
int maxx = 0;
for (int i = 1; i <= n; i++)
	for (int j = 1; j <= n; j++) {
		int x_1 = i - h[i][j] + 1, y_1 = lmin[i][j] + 1;
		int x_2 = i, y_2 = rmin[i][j] - 1;
		int t = cut(x_1, y_1, x_2, y_2);
		if (t > maxx) {
			maxx = t;
			s.first = x_1;
			s.second = y_1;
			e.first = x_2;
			e.second = y_2;
		}
	}

5.一些判定

if (maxx >= k && maxx <= 2 * k) {
	cout << s.second << " " << s.first << " " << e.second << " " << e.first;
	return 0;
}
if (maxx < k) {
	cout << "NIE";
	return 0;
}

6.逐行判定并分割极大子矩阵来求得答案（这几个if也可以不写因为前两个里都含有return 0，我写上是为了更直观）

for (int i = s.first; i <= e.first; i++) {
	int t = cut(i, s.second, i, e.second);
	if (t >= k && t <= 2 * k) {
		cout << s.second << " " << i << " " << e.second << " " << i;
		return 0;
	}
	if (t > 2 * k) {
		int sum = 0;
		for (int j = s.second; j <= e.second; j++) {
			if (sum + a[i][j] > 2 * k) {
				cout << s.second << " " << i << " " << j - 1 << " " << i;
				return 0;
			}
			sum += a[i][j];
		}
	}
	if (t < k) {
		maxx -= t;
		if (maxx >= k && maxx <= 2 * k) {
			cout << s.second << " " << i + 1 << " " << e.second << " " << e.first;
			return 0;
		}
	}
}

完整代码：

注释就不写了，之前已经讲得很详细了

#include <bits/stdc++.h>
using namespace std;
int n, a[2001][2001], lmin[2001][2001], rmin[2001][2001], h[2001][2001];
long long cnt[2001][2001];
void solve(int line) {
	stack<int> l, r;
	for (int i = 1; i <= n; i++) {
		while (!l.empty() && h[line][l.top()] >= h[line][i])
			l.pop();
		if (l.empty())
			lmin[line][i] = 0;
		else
			lmin[line][i] = l.top();
		l.push(i);
	}
	for (int i = n; i >= 1; i--) {
		while (!r.empty() && h[line][r.top()] >= h[line][i])
			r.pop();
		if (r.empty())
			rmin[line][i] = n + 1;
		else
			rmin[line][i] = r.top();
		r.push(i);
	}
}
long long cut(int x_1, int y_1, int x_2, int y_2) {
	return cnt[x_2][y_2] - cnt[x_2][y_1 - 1] - cnt[x_1 - 1][y_2] + cnt[x_1 - 1][y_1 - 1];
}
int main() {
	int k;
	cin >> k >> n;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++) {
			cin >> a[i][j];
			if (a[i][j] >= k && a[i][j] <= 2 * k) {
				cout << j << " " << i << " " << j << " " << i;
				return 0;
			}
			cnt[i][j] = cnt[i - 1][j] + cnt[i][j - 1] - cnt[i - 1][j - 1] + a[i][j];
		}
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++)
			if (a[i][j] <= 2 * k)
				h[i][j] = h[i - 1][j] + 1;
	for (int i = 1; i <= n; i++)
		solve(i);
	pair<int, int> s, e;
	int maxx = 0;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++) {
			int x_1 = i - h[i][j] + 1, y_1 = lmin[i][j] + 1;
			int x_2 = i, y_2 = rmin[i][j] - 1;
			int t = cut(x_1, y_1, x_2, y_2);
			if (t > maxx) {
				maxx = t;
				s.first = x_1;
				s.second = y_1;
				e.first = x_2;
				e.second = y_2;
			}
		}
	if (maxx >= k && maxx <= 2 * k) {
		cout << s.second << " " << s.first << " " << e.second << " " << e.first;
		return 0;
	}
	if (maxx < k) {
		cout << "NIE";
		return 0;
	}
	for (int i = s.first; i <= e.first; i++) {
		int t = cut(i, s.second, i, e.second);
		if (t >= k && t <= 2 * k) {
			cout << s.second << " " << i << " " << e.second << " " << i;
			return 0;
		}
		if (t > 2 * k) {
			int sum = 0;
			for (int j = s.second; j <= e.second; j++) {
				if (sum + a[i][j] > 2 * k) {
					cout << s.second << " " << i << " " << j - 1 << " " << i;
					return 0;
				}
				sum += a[i][j];
			}
		}
		if (t < k) {
			maxx -= t;
			if (maxx >= k && maxx <= 2 * k) {
				cout << s.second << " " << i + 1 << " " << e.second << " " << e.first;
				return 0;
			}
		}
	}
	return 0;
}