【BZOJ1478】Sgu282 Isomorphism Pólya定理神题

【BZOJ1478】Sgu282 Isomorphism

题意：用$m$种颜色去染一张$n$个点的完全图，如果一个图可以通过节点重新标号变成另外一个图，则称这两个图是相同的。问不同的染色方案数。答案对$P$取模。
$nle 53,mle 1000,P>n,P$是质数。
题解：对于本题来说，每个元素是所有边，每个置换是边的置换，而边的置换难以表示，点的置换容易表示，所以我们考虑点置换和边置换的关系。
如果两个点置换有着相同的结构，则它们对应的边置换的循环数相同。
$$egin{pmatrix}1&2&3 2&1&3end{pmatrix} ightarrowegin{pmatrix}(1,2)&(1,3)&(2,3) (1,2)&(2,3)&(1,3)end{pmatrix}$$
$$egin{pmatrix}1&2&3 1&3&2end{pmatrix} ightarrowegin{pmatrix}(1,2)&(1,3)&(2,3) (1,3)&(1,2)&(2,3)end{pmatrix}$$
一个点置换的结构：在表示成循环节时，循环的大小依次为$L_1,L_2,...L_K$。
$L_1,L_2,...L_K$可以组成一个$n$的划分。
而当$n$=53时，$n$个划分数不超过$300000$。
可以通过搜索得出。
那么一个结构为$L_1,L_2,...L_K$的点置换对应边置换的循环数是什么呢？
对于端点在不同点循环中的边，设两端点所在点循环大小为$L_1,L_2$，则这样的边一共有$L_1 imes L_2$条，而每个循环有$lcm(L_1,L_2)$条边，所以一共有$gcd(L_1,L_2)$个循环。
对于端点在同一点循环中的边，设所在点循环的大小为$L$，这样的边一共有$C_L^2$条。然后分奇偶讨论：

1. 如果$L$是奇数，那么一个循环中有$L$条边，所以循环数为${C_L^2over L}={L-1over 2}$。
2. 如果$L$是偶数，那么一个循环中有$L$条边，但是如果两端点相隔$Lover 2$，这个循环中有$Lover 2$条边。所以循环数为${C_L^2-{Lover 2}over L}+1={Lover 2}$。

所以一个大小为$L$的点循环中边循环的数量是$lfloor{Lover 2} floor$。
结构为$L_1,L_2,...L_K$的点置换中边循环的数量就是
$$sumlimits_{i=1}^nlfloor{L_iover 2} floor+sumlimits_{i=1}^{Ksumlimits_{j=1}}{i-1}gcd(L_i,L_j)$$
那么有多少结构为$L_1,L_2...L_K$的点置换呢？可以先求出$n!over {L_1!L_2!...L_K!}$得到每个点属于哪个循环的方案数，而对于每个循环，我们可以先确定第一个点，然后其余点任意排列，方案数是$(L_1-1)!(L_2-1)!...(L_K-1)!$。乘起来得到$n!over {L_1L_2...L_k}$。
又由于存在一些$L$相等的情况，设有$t$种本质不同的$L$，每种个数为$B_i$，所以总方案数还要除掉$B_1!B_2!...B_t!$，所以总数为：
$$n!over L_1L_2...L_KB_1!B_2!...B_t!$$
最后套上Pólya定理，令$c=sumlimits_{i=1}^nlfloor{L_iover 2} floor+sumlimits_{i=1}^{Ksumlimits_{j=1}}{i-1}gcd(L_i,L_j)$，$s={n!over L_1L_2...L_KB_1!B_2!...B_t!}$，答案就是
$$frac 1 {n!}sum s imes m^c$$

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
typedef long long ll;
ll P,ans;
int n,m;
ll jc[70],jcc[70],ine[70],pw[3700];
int g[70][70],v[70];
inline int rd()
{
	int ret=0,f=1;	char gc=getchar();
	while(gc<'0'||gc>'9')	{if(gc=='-')	f=-f;	gc=getchar();}
	while(gc>='0'&&gc<='9')	ret=ret*10+(gc^'0'),gc=getchar();
	return ret*f;
}
inline ll C(int a,int b)
{
	if(a<b)	return 0;
	return jc[a]*jcc[a-b]%P*jcc[b]%P;
}
void dfs(int len,int sum)
{
	if(!sum)
	{
		ll s=jc[n],c=0;
		int i,j,tmp;
		for(i=1,tmp=0;i<=len;i++)
		{
			tmp++;
			if(i==len||v[i]!=v[i+1])	s=s*jcc[tmp]%P,tmp=0;
			s=s*ine[v[i]]%P;
			c+=v[i]>>1;
			for(j=1;j<i;j++)	c+=g[v[j]][v[i]];
		}
		ans=(ans+pw[c]*s)%P;
		return ;
	}
	for(int i=min(sum,v[len]);i>=1;i--)	v[len+1]=i,dfs(len+1,sum-i);
}
int main()
{
	n=rd(),m=rd(),P=rd();
	int i,j;
	for(pw[0]=i=1;i<=n*n/2;i++)	pw[i]=pw[i-1]*m%P;
	for(i=0;i<=n;i++)	g[i][0]=i;
	for(i=0;i<=n;i++)	for(j=1;j<=i;j++)	g[i][j]=g[j][i%j];
	jc[1]=ine[1]=jcc[1]=jc[0]=ine[0]=jcc[0]=1;
	for(i=2;i<=n;i++)	jc[i]=jc[i-1]*i%P,ine[i]=P-(P/i)*ine[P%i]%P,jcc[i]=jcc[i-1]*ine[i]%P;
	v[0]=n;
	dfs(0,n);
	printf("%lld",ans*jcc[n]%P);
	return 0;
}