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

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

---

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

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

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

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

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

By the way，听说此题有人直接算$\frac{C_{Sr+Sb}^{Sr}*C_{Sr+Sb+Sg}^{Sg}}{m+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;
}