好题！感觉比那些写出 DP 然后无脑上 GF 数学方法硬推的题有价值。
首先有个朴素的想法：设 \(f_{l, r, u, d}\) 表示这个矩形的方案数，那么枚举分界点转移。
引用大佬的话：
直接枚举分界点显然会算重，这里考虑有以下两种一般性的思路：
- 探究一下什么序列能够被得到，由此而不是由操作序列的角度设计 DP 状态。
- 硬是从操作序列的角度入手，不过给操作序列定个顺序（比方说，字典序最小的操作序列），以保证任意一个合法的序列只会被算一次。
发现能得到的序列没啥规律，我们考虑给操作序列定顺序。我们将对 \(x=i\) 进行的操作记作 \(x(i)\)，把对 \(y=i\) 进行的操作记作 \(y(i)\)，那么我们定义操作的大小关系为 \(x(1) < y(1) < x(2) < y(2) < \cdots < x(k) < y(k)\)（\(k\) 为当前矩形内的点数）。
这样我们重新定义 DP：设 \(f_{l, r, u, d}\) 表示这个矩形内由字典序最小的操作序列得到的方案数，这样能够保证方案不会算重，因为只有字典序最小的被计算了一次。
字典序最小我们肯定贪心的使得操作第一步最小，所以我们先枚举第一步操作。
发现我们没法很好的保证字典序最小，那么可以考虑容斥。我们可以将第一步操作不是最小，但是能够拼出当前序列的方案数减去。
假如我们现在考虑第一步操作为 \(x(i)\)，那么假如左边的矩形内第一步操作为 \(x(j)\)，那么我们需要保证先 \(x(j)\) 后 \(x(i)\) 不能得到当前的序列，那么就减去先 \([l, j - 1]\) 再 \([j, i - 1]\) 最后 \([i, r]\) 能得到的字典序最小的方案数。
我们可以记录一个 \(fx_i\) 为第一步操作为 \(x(i)\) 时，字典序最小的左边矩形的方案数，这样上述容斥就是 \(fx_i = f_{l, i, u, d} - \sum_{j < i}fx_j \times f_{j, i - 1, u, d}\)。
再考虑左边矩形第一步操作为 \(y(j)\) 的情况，那么如果 \(y(j)\) 能先执行，说明 \(y(j),x(i)\) 划分的四个矩形中右上角的矩形为空（否则 \(y(j)\) 划分必定先会右上角，不可能与先划分 \(x(i)\) 相同），那么我们就只考虑右上角为空的 \(y(j)\)，然后类似的设 \(fy_i\) 进行转移，\(fx_i \gets fx_i - \sum_{j \le \min\{p_i,\cdots,p_r\}}fy_j \times f_{l, i - 1, j, d}\)。
那么 \(fx_i\) 对答案造成的贡献就是 \(f_{l, r, u, d} \gets f_{l, r, u, d} + fx_i \times f_{i, r, u, d}\)。
\(fy_i\) 进行类似的转移即可。
这样状态数为 \(O(n^4)\)，转移为 \(O(n^2)\)，总复杂度为 \(O(n^6)\)。
同类型的题可以看 AGC056B Range Argmax，也是通过确定操作序列的顺序使得操作序列与结果一一对应。
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 55, P = 1000000007;
int n, p[MAXN], q[MAXN];
int f[MAXN][MAXN][MAXN][MAXN];
bool inc(int x, int l, int r) {
return l <= x && x <= r;
}
int dp(int l, int r, int u, int d) {
vector<int> x(MAXN), y(MAXN), xs(MAXN), ys(MAXN), pi(MAXN), qi(MAXN), pim(MAXN), qim(MAXN);
if (l > r || u > d) return 0;
if (f[l][r][u][d] != -1) return f[l][r][u][d];
int res = 0;
f[l][r][u][d] = 0;
int m = 0;
for (int i = l; i <= r; i++)
if (inc(p[i], u, d)) ++m, xs[m] = i, ys[m] = p[i];
sort(ys.begin() + 1, ys.begin() + 1 + m);
for (int i = 1; i <= m; i++) {
pi[i] = lower_bound(ys.begin() + 1, ys.begin() + 1 + m, p[xs[i]]) - ys.begin();
qi[pi[i]] = i;
}
if (m <= 1) return f[l][r][u][d] = 1;
if (!inc(p[l], u, d)) return f[l][r][u][d] = dp(l + 1, r, u, d);
if (!inc(p[r], u, d)) return f[l][r][u][d] = dp(l, r - 1, u, d);
if (!inc(q[u], l, r)) return f[l][r][u][d] = dp(l, r, u + 1, d);
if (!inc(q[d], l, r)) return f[l][r][u][d] = dp(l, r, u, d - 1);
pim[m] = pi[m], qim[m] = qi[m];
for (int i = m - 1; i >= 1; i--) {
pim[i] = min(pim[i + 1], pi[i]);
qim[i] = min(qim[i + 1], qi[i]);
}
for (int i = 1; i <= m; i++) {
x[i] = dp(l, xs[i] - 1, u, d);
for (int j = 1; j < i; j++) {
x[i] = (x[i] - 1ll * x[j] * dp(xs[j], xs[i] - 1, u, d) % P + P) % P;
}
for (int j = 1; j <= pim[i]; j++) {
x[i] = (x[i] - 1ll * y[j] * dp(l, xs[i] - 1, ys[j], d) % P + P) % P;
}
res = (res + 1ll * x[i] * dp(xs[i], r, u, d)) % P;
y[i] = dp(l, r, u, ys[i] - 1);
for (int j = 1; j < i; j++) {
y[i] = (y[i] - 1ll * y[j] * dp(l, r, ys[j], ys[i] - 1) % P + P) % P;
}
for (int j = 1; j <= qim[i]; j++) {
y[i] = (y[i] - 1ll * x[j] * dp(xs[j], r, u, ys[i] - 1) % P + P) % P;
}
res = (res + 1ll * y[i] * dp(l, r, ys[i], d)) % P;
}
return f[l][r][u][d] = res;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
scanf("%d", &p[i]);
q[p[i]] = i;
}
memset(f, -1, sizeof f);
printf("%d\n", dp(1, n, 1, n));
return 0;
}