题目链接：BZOJ - 1004
题目分析
首先，几个定义和定理引理：
群：G是一个集合，*是定义在这个集合上的一个运算。
如果满足以下性质，那么(G, *)是一个群。
1）封闭性，对于任意 a, b 属于 G， a * b 属于 G
2）结合律， a * b * c = a * (b * c)
3）单位元，在 G 中存在一个单位元 e ，使得对于 G 中任意的 a ， a * e = e * a = a
4）逆元， 对于 G 中任意的 a ，在 G 中存在 b ， 使得 a * b = e , 其中 b 叫做 a 的逆元
比如在模一个数意义下的整数加法就是一个群。
满足交换律的群是交换群，又叫阿贝尔群。
置换：可以用 (a1 -> b1, a2 -> b2, ... , an -> bn) 表示一个置换，其中 a1, ... , an 和 b1, ..., bn 都是1 到 n 的一个排列；
如果一些置换和它们的叠加运算构成一个群，就把它们叫做一个置换群。
在置换群中的 Burnside 引理：如果按照一定要求，要对1到n 的位置染色，那么本质不同的染色方案数为置换群中每个置换的不动染色方案数的平均数。
来解释以下，本质不同的染色方案是指，两个染色方案不能通过置换群中的任意置换变换使其相同，那么它们就是本质不同的。
某个置换的不动染色方案数是指，用这个置换变换之后没有发生变化的染色方案。
那么我们就是要求出每个置换的不动染色方案数。
Polya定理告诉我们，如果是用k种颜色染色，那么对于置换 P 来说，它的不动染色方案数为 k^(L(P)), 其中L(P)为置换P的循环节数。
这个很好理解，因为每个循环节如果要使置换之后不改变的话，应该使这个循环节中都染同一种颜色。
但是，这道题的染色要求是给定了3种颜色的数量，就不能用 k^(L(P)) 来求了。
我们使用一种3维DP，还是先搜索一遍求出这个置换的所有循环节的大小。然后把每个循环节看作一个物品，染一种颜色。
用 f[i][j][k] 表示染了前 i 件物品，第1种颜色用了 j 个， 第2种颜色用了 k 次的染色方案数。
这道题还需要注意的一点是，需要加上一个不做任何改变的置换作为单位元才能成为一个置换群。
这个单位元的不动染色方案数就是所有的染色方案数，直接用多重集排列就好了， n! / (a! b! c!)。
代码
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
const int MaxN = 60;
int Na, Nb, Nc, n, m, Mod, Ans, Top;
int A[MaxN], V[MaxN], Sum[MaxN], f[MaxN][MaxN][MaxN];
bool Used[MaxN];
inline int Pow(int a, int b)
{
int ret = 1, f = a;
while (b)
{
if (b & 1)
{
ret *= f;
ret %= Mod;
}
b >>= 1;
f *= f;
f %= Mod;
}
return ret;
}
inline int NY(int x) {return Pow(x, Mod - 2);}
void DFS(int x, int Num)
{
if (Used[x])
{
V[++Top] = Num;
return;
}
Used[x] = true;
DFS(A[x], Num + 1);
}
inline int Fac(int x)
{
int ret = 1;
for (int i = 1; i <= x; ++i)
{
ret *= i;
ret %= Mod;
}
return ret;
}
int main()
{
scanf("%d%d%d%d%d", &Na, &Nb, &Nc, &m, &Mod);
n = Na + Nb + Nc;
Ans = Fac(n) * NY(Fac(Na)) % Mod * NY(Fac(Nb)) % Mod * NY(Fac(Nc)) % Mod;
for (int i = 1; i <= m; ++i)
{
for (int j = 1; j <= n; ++j)
{
scanf("%d", &A[j]);
Used[j] = false;
}
Top = 0;
for (int j = 1; j <= n; ++j)
if (!Used[j]) DFS(j, 0);
Sum[0] = 0;
for (int j = 1; j <= Top; ++j)
Sum[j] = Sum[j - 1] + V[j];
f[0][0][0] = 1;
for (int j = 1; j <= Top; ++j)
{
for (int p = 0; p <= Na; ++p)
{
for (int q = 0; q <= Nb; ++q)
{
f[j][p][q] = 0;
if (p >= V[j])
{
f[j][p][q] += f[j - 1][p - V[j]][q];
f[j][p][q] %= Mod;
}
if (q >= V[j])
{
f[j][p][q] += f[j - 1][p][q - V[j]];
f[j][p][q] %= Mod;
}
if ((Sum[j] - p - q) >= V[j])
{
f[j][p][q] += f[j - 1][p][q];
f[j][p][q] %= Mod;
}
}
}
}
Ans += f[Top][Na][Nb];
Ans %= Mod;
}
Ans *= NY(m + 1);
Ans %= Mod;
printf("%d
", Ans);
return 0;
}