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

题目传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=1004

题目分析：最近我总是在省选题中找题目做，结果遇到各种奇奇怪怪的坑。

这题的题面告诉我们，如果某个状态可以通过一次洗牌到达另一种状态，就从这种状态向另一种状态连一条有向边，那么最终的图一定全部是双向边，并且所有的状态组成了一些团。团的个数即为答案。现在的问题是如何求出团的个数？

我自己的想法是用2m枚举一些洗牌，然后看看有多少个状态通过这些洗牌还是洗回自身（这一点可以通过并查集+带权01背包知道），然后跑个容斥，就可以知道向外连出x条边的状态有多少个，记为num[x]。由于所有状态组成一些团，它对答案的贡献即为num[x]x+1。

然而这样显然超时，我想了很久，然后看了题解才知道这里要用到一个叫做Burnside引理的东西，它貌似和Polya定理有些关联。于是我果断%了一波网上dalao写的关于这个引理的学习笔记。虽然证明看得差不多了，但其实还有两个核心的地方不是很明白：同一个等价类的元素，Zi肯定是相同的，以及i的等价类的个数乘以令i不变化的置换的个数=置换的总个数。在此留一个坑，希望日后自己能来填或者有大神路过教导一下QAQ。

上述引理告诉我们：对于一个置换群，只要求出每一个置换的不动点个数，其平均值就是该群的本质不同的染色方案数。也就是说在这题中，我们只要对每一个洗牌方案算出有多少种状态，洗了之后还是它自身，加起来除以m+1即可。为什么是m+1呢？因为我们还要算上这个群的单位元——p[i]=i的置换，它对答案的贡献为CSrSr+Sb∗CSgSr+Sb+Sg（但其实直接用DP来算即可）。除以m+1时算个逆元。

By the way，听说此题有人直接算CSrSr+Sb∗CSgSr+Sb+Sgm+1就过了，但我对这个有些不解：每个团的大小一定是m+1吗原谅我不会证？难道题目有什么特殊保证而我没有看出来？总之本人是写了DP的。

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=25;
const int maxm=65;

int f[maxn][maxn][maxn];
int X,Y;

int Size[maxm];
int fa[maxm];

int g[maxm][maxm];
int a,b,c,m,p;
int n;

int Up(int x)
{
if (fa[x]==x) return x;
return (fa[x]=Up(fa[x]));
}

void Add(int x,int y)
{
x=Up(x);
y=Up(y);
if (x==y) return;
fa[x]=y;
Size[y]+=Size[x];
}

void Exgcd(int a,int b)
{
if (!b) X=1,Y=0;
else
{
Exgcd(b,a%b);
int u=Y;
int v=X-a/b*Y;
X=u;
Y=v;
}
}

int main()
{
freopen("1004.in","r",stdin);
freopen("1004.out","w",stdout);

scanf("%d%d%d%d%d",&a,&b,&c,&m,&p);
n=a+b+c;
for (int i=1; i<=m; i++)
for (int j=1; j<=n; j++) scanf("%d",&g[i][j]);
m++;
for (int i=1; i<=n; i++) g[m][i]=i;

int ans=0;
for (int i=1; i<=m; i++)
{
for (int j=1; j<=n; j++) fa[j]=j,Size[j]=1;
for (int j=1; j<=n; j++) Add(j,g[i][j]);
for (int j=1; j<=n; j++) Up(j);

memset(f,0,sizeof(f));
f[0][0][0]=1;
for (int j=1; j<=n; j++)
if (fa[j]==j)
{
int k=Size[j];
for (int x=a; x>=0; x--)
for (int y=b; y>=0; y--)
for (int z=c; z>=0; z--)
{
int &v=f[x][y][z];
if (x>=k) v=(v+f[x-k][y][z])%p;
if (y>=k) v=(v+f[x][y-k][z])%p;
if (z>=k) v=(v+f[x][y][z-k])%p;
}
}
ans=(ans+f[a][b][c])%p;
}

Exgcd(p,m);
while (Y<0) Y+=p;
ans=(ans*Y)%p;
printf("%d\n",ans);

return 0;
}