[BZOJ1004][HNOI2008]Cards（Burnside引理+DP）

加上一个置换$1\>2\>3\>...\>n$，形成$m+1$个置换。
根据Burnside引理，答案等于：

$$\frac{1}{m+1}\sum_{x=1}^{m+1}D(a_x)$$
$D(a_x)$表示在置换 $a_x$（第 $x$个置换）下，不变的元素（染色方案）个数。
对于每个$1\leq x\leq m+1$，怎么求 $D(a_x)$呢？
首先把 $a_x$分循环节。可以得出，一种方案在置换 $a_x$下不变，等价于这种方案在 $a_x$的每个循环节内只有一种颜色。
如一种置换为 $2\>4\>5\>1\>3$，循环节为 $(1,2,4)(3,5)$，那么一种方案在置换 $a_x$下不变，等价于这种方案的第 $1$张，第$2$张和第 $4$张牌同色，第$3$张牌和第 $5$张牌同色。
然后？DP！
设$f[k][i][j]$表示到了第 $k$个循环节，用了$i$张红色 $j$张蓝色的方案数。边界$f[0][0][0]=1$。转移就是枚举第 $k$个循环节被染成的颜色（$tot_k$表示第 $k$个循环节包含的牌数）：
$$f[k][i][j]=f[k-1][i-tot_k][j]+f[k-1][i][j-tot_k]+f[k-1][i][j]$$
这样， $D(a_x)$就等于 $f[a_x$的循环节个数 $][Sr][Sb]$。
代码：

#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 N = 66, M = 23;
int n, Sr, Sb, Sg, m, MX, to[N][N], f[N][M][M];
bool vis[N];
int qpow(int a, int b) {
    int res = 1; while (b) b & 1 ? res = res * a % MX : 0,
        a = a * a % MX, b >>= 1; return res;
}
int main() {
    int i, j, k, h; n = (Sr = read()) + (Sb = read()) + (Sg = read());
    m = read(); MX = read(); for (i = 1; i <= m; i++) for (j = 1; j <= n; j++)
        to[i][j] = read(); m++; for (j = 1; j <= n; j++) to[m][j] = j;
    int ans = 0; for (i = 1; i <= m; i++) {
        memset(vis, 0, sizeof(vis)); memset(f, 0, sizeof(f)); f[0][0][0] = 1;
        int t = 0; for (j = 1; j <= n; j++) if (!vis[j]) {
            int x = to[i][j], tot = 1; vis[j] = 1; t++;
            while (x != j) vis[x] = 1, x = to[i][x], tot++;
            for (k = 0; k <= Sr; k++) for (h = 0; h <= Sb; h++) {
                if (k >= tot)
                    f[t][k][h] = (f[t][k][h] + f[t - 1][k - tot][h]) % MX;
                if (h >= tot)
                    f[t][k][h] = (f[t][k][h] + f[t - 1][k][h - tot]) % MX;
                f[t][k][h] = (f[t][k][h] + f[t - 1][k][h]) % MX;
            }
        }
        ans = (ans + f[t][Sr][Sb]) % MX;
    }
    cout << ans * qpow(m, MX - 2) % MX << endl;
    return 0;
}