[BZOJ5010][Fjoi2017]矩阵填数（状压DP）

从最简单的问题开始考虑：
一个面积为$s$的矩阵，每个位置可以填上$[1,m]$中的任意一个整数，并且要求这个矩阵的最大值为$v$，求方案数。
这个问题等价于每个位置填上$[1,v]$中的任意一个整数，并且至少要有一个$v$。
很容易得出答案为$v^s-(v-1)^s$。
而如果这$n$个矩阵两两没有交集，那么可以简单地把这个矩阵分成$n+1$个部分，分别计算结果之后相乘。
下面就样例中的第$1$组数据，讨论有交集的情况。
可以看出，第$1$组数据中的$2$个矩形把矩阵分成了$4$个部分：
1、$(1,3)(3,1)$
2、$(1,1)(1,2)(2,1)$
3、$(2,3)(3,2)(3,3)$
4、$(2,2)$
第$1$个部分无限制，每个格子可以填$1$或$2$，方案数为$4$。
第$4$个部分恰好是两个矩形的交集，因此绝对不能填大于$1$的数，否则第$2$个矩形的最大值永远不会等于$1$。因此$(2,2)$只能填$1$。
第$3$个部分只能填$1$，方案数为$1$。
第$2$个部分可以填$1$或$2$，且要至少有一个$2$，所以方案数为$2^3-(2-1)^3=7$。
综上所述，方案数为$4\times 1\times 1\times 7=28$。
而在实际上，$n$个矩形可以把矩阵分割成$O(2^n)$块。
于是就有了一个状压DP的模型：$f[i,S]$表示到了第$i$块，已经满足限制的子矩形集合为$S$的方案数。当然，这些块不包括子矩形外面无限制的部分（最后乘上无限制部分的方案数即可）。
首先需要预处理出每块是哪些子矩形的交集，每一块的面积$area_i$，以及这些子矩形需求的最大值的最小值$val$。也就是第$i$块的最大值不能超过$val_i$。
还要对于每块预处理出一个集合$cov_i$，表示包含这一块的所有子矩形中，需求的最大值等于$val_i$的子矩形集合。
DP转移：
1、填小于$val_i$的数：

$$f[i+1][S]+=f[i][j]\times(val_{i+1}-1)^{area_{i+1}}$$
2、填 $[1,val_i]$的数，且至少一个 $val_i$：
$$f[i+1][S|cov_{i+1}]+=f[i][j]\times(val_{i+1}^{area_{i+1}}-(val_{i+1}-1)^{area_{i+1}})$$
理论复杂度 $O(T\times 4^n)$，但实际上远不达这个上界。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int E = 13, N = 1100, MX = 1e9 + 7;
int h, w, m, n, f[N][N], Cm, ara[N], cov[N], tot, sta[N];
int qpow(int a, int b) {
    int res = 1; while (b) b & 1 ? res = 1ll * res * a % MX : 0,
        a = 1ll * a * a % MX, b >>= 1; return res;
}
struct cyx {
    int X1, Y1, X2, Y2, V; cyx() {}
    cyx(int _, int __, int ___, int ____, int _____) :
        X1(_), Y1(__), X2(___), Y2(____), V(_____) {}
    int area() {return (X2 - X1 + 1) * (Y2 - Y1 + 1);}
} orz[E], uni[N];
cyx mers(cyx a, cyx b) {
    int l1 = a.X1, r1 = a.X2, l2 = b.X1, r2 = b.X2; cyx res;
    if (l1 > l2) swap(l1, l2), swap(r1, r2); if (r1 < l2) return
    cyx(0, 0, 0, 0, 0); res.X1 = l2; res.X2 = min(r1, r2);
    l1 = a.Y1; r1 = a.Y2; l2 = b.Y1; r2 = b.Y2;
    if (l1 > l2) swap(l1, l2), swap(r1, r2); if (r1 < l2) return
    cyx(0, 0, 0, 0, 0); res.Y1 = l2; res.Y2 = min(r1, r2);
    res.V = min(a.V, b.V); return res;
}
void dfs(int dep, cyx s, int st) {
    if (dep == n + 1) return (void) (uni[st] = s);
    dfs(dep + 1, s, st); cyx t = mers(s, orz[dep]);
    if (t.X1) dfs(dep + 1, t, st | (1 << dep - 1));
}
void work() {
    int i, j; h = read(); w = read(); m = read(); n = read();
    for (i = 1; i <= n; i++) orz[i].X1 = read(), orz[i].Y1 = read(),
        orz[i].X2 = read(), orz[i].Y2 = read(), orz[i].V = read();
    Cm = 1 << n; for (i = 0; i < Cm; i++) uni[i] = cyx(1, 1, 0, 0, m);
    dfs(1, cyx(1, 1, h, w, m), 0); for (i = 0; i < Cm; i++)
        ara[i] = uni[i].area(); for (i = Cm - 1; i; i--) {
        ara[0] -= ara[i]; for (j = (i - 1) & i; j; j = (j - 1) & i)
            ara[j] -= ara[i];
    }
    for (i = 0; i < Cm; i++) cov[i] = 0;
    for (i = 0; i < Cm; i++) for (j = 1; j <= n; j++)
        if (((i >> j - 1) & 1) && orz[j].V == uni[i].V)
            cov[i] |= 1 << j - 1; tot = 0;
    for (i = 1; i < Cm; i++) if (ara[i]) sta[++tot] = i;
    for (i = 0; i <= tot; i++) for (j = 0; j < Cm; j++) f[i][j] = 0; f[0][0] = 1;
    for (i = 0; i < tot; i++) {
        int S = sta[i + 1], xs = qpow(uni[S].V - 1, ara[S]),
            ys = (qpow(uni[S].V, ara[S]) - xs + MX) % MX;
        for (j = 0; j < Cm; j++) if (f[i][j])
            f[i + 1][j] = (f[i + 1][j] + 1ll * f[i][j] * xs % MX) % MX,
            f[i + 1][j | cov[S]] = (f[i + 1][j | cov[S]] + 1ll *
                f[i][j] * ys % MX) % MX;
    }
    printf("%d\n", 1ll * f[tot][Cm - 1] * qpow(m, ara[0]) % MX);
}
int main() {
    int T = read(); while (T--) work();
    return 0;
}