传送门：http://www.lydsy.com:808/JudgeOnline/problem.php?id=3294

这道题的dp方程容易想到：令f[i][j][k]表示前i种颜色占了j行，k列的方案数，g[i][j][k] 表示用第i种颜色占了j行，k列的方案总数，则f[i][j][k] = sigma(f[i - 1][x][y] * g[i][j - x][k - y])；

关键是g[i][j][k]怎么求，这就是经典的容斥原理了。

先减去少一行或一列的，再加上少两行或少两列或少一行一列的......

 

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long LL;
const LL mod = 1000000009LL;
const int maxn = 35, maxk = 15;
LL f[maxk][maxn][maxn];
LL g[maxk][maxn][maxn];
LL c[maxn * maxn][maxn * maxn];
int d[maxk];
int n, m, k;
LL ans;
void pre() {
    c[0][0] = 1;
    for(int i = 1; i < maxn * maxn; i ++) {
        c[i][0] = c[i][i] = 1;
    }
    for(int i = 2; i < maxn * maxn; i ++) {
        for(int j = 1; j < i; j ++) {
            c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
        }
    }
    return;
}
void cal() {
    for(int i = 1; i <= k; i ++) {
        for(int a = 1; a <= n; a ++) {
            for(int b = 1; b <= m; b ++) {
                for(int x = 1; x <= a; x ++) {
                    for(int y = 1; y <= b; y ++) {
                        if(((a + b) & 1) == ((x + y) & 1)) {
                            g[i][a][b] += c[x * y][d[i]] * c[a][x] % mod * c[b][y] % mod;
                        } else {
                            g[i][a][b] -= c[x * y][d[i]] * c[a][x] % mod * c[b][y] % mod;
                        }
                        if(g[i][a][b] < 0) g[i][a][b] += mod;
                        g[i][a][b] %= mod;
                    }
                }
            }
        }
    }
    return;
}
int main() {
    pre();
    scanf("%d%d%d", &n, &m, &k);
    for(int i = 1; i <= k; i ++) {
        scanf("%d", &d[i]);
    }
    cal();
    f[0][0][0] = 1;
    for(int i = 1; i <= k; i ++) {
        for(int a = i; a <= n; a ++) {
            for(int b = i; b <= m; b ++) {
                for(int x = max(a - d[i], 0); x < a; x ++) {
                    for(int y = max(b - d[i], 0); y < b; y ++) {
                        f[i][a][b] = (f[i][a][b] + f[i - 1][x][y] * g[i][a - x][b - y] % mod * c[n - x][a - x] % mod * c[m - y][b - y] % mod) % mod;
                    }
                }
            }
        }
    }
    for(int i = 0; i <= n; i ++) {
        for(int j = 0; j <= m; j ++) {
            ans = (ans + f[k][i][j]) % mod;
        }
    }
    cout << ans << endl;
    return 0;
}