【BZOJ1488】[HNOI2009] 图的同构（无向图Polya定理）

点此看题面

大致题意： 求含\(n\)个点的互不同构的简单无向图个数。

转化

考虑把一条边是否存在看成黑白染色，这就变成了一个经典的\(Polya\)计数问题。

应该可以看作是无向图\(Polya\)定理的板子吧，但确实挺难的。

无向图\(Polya\)定理

本题中置换群的对象就是这\(\frac{n(n-1)}2\)条边，但由于置换的数量达到\(n!\)，显然不可以裸暴力。

假设点的置换是\(i\rightarrow P_i\)，则边的置换就是\((i,j)\rightarrow(P_i,P_j)\)。

考虑两种置换循环节个数的关系，发现：

• 若点\(i,j\)属于同一长度为\(x\)的循环中，\((i,j)\)组成的置换中循环节个数为\(\lfloor\frac x2\rfloor\)。
• 若点\(i,j\)分别属于长度为\(x,y\)的两个循环中，\((i,j)\)组成的置换中循环节个数为\(gcd(x,y)\)。

因此，如果一个点置换长度分别为\(L_1,L_2,...,L_t\)，则边置换的循环节总数就是：

\[\sum_{i=1}^t\lfloor\frac{L_i}2\rfloor+\sum_{i=1}^t\sum_{j=i+1}^tgcd(L_i,L_j) \]

那么就需要求满足循环节长度为\(L_1,L_2,...,L_t\)的边置换数目，也就是把\(1\sim n\)放入大小分别为\(L_1,L_2,...,L_t\)的环的方案数。

如果\(L\)互不相同，方案数就是：

\[\frac{n!}{\prod_{i=1}^tL_i} \]

而若可能有相同的\(L\)，同样大小的环就会造成重复贡献。

因此我们将\(L_i\)去重，用\(C_i\)表示第\(i\)种值的个数，方案数就是：

\[S=\frac{n!}{\prod_{i=1}^t(L_i)^{C_i}\times C_i!} \]

在这种表示方法下，之前的循环节总数就应该是：

\[T=\sum_{i=1}^tC_i\times \lfloor\frac{L_i}2\rfloor+\sum_{i=1}^t(\frac{C_i\times (C_i-1)}2\times L_i+\sum_{j=t+1}^n C_i\times C_j\times gcd(L_i,L_j)) \]

最终答案就是：（注意这其实是一个\(Polya\)定理的形式，只是在\(m^{c(g_i)}\)之前乘上了一个方案数）

\[ans=\frac1{n!}\sum S\times 2^T \]

这道题的核心思想就是把相似的情况放在一起讨论，除去冗余的方法在\(Polya\)计数问题中是一个很有用的优化。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 60
#define X 997
using namespace std;
int n,ans,t,l[N+5],c[N+5],pw[X],Inv[N+5],Fac[N+5],IFac[N+5],g[N+5][N+5];
I int QP(RI x,RI y) {RI t=1;W(y) y&1&&((t*=x)%=X),(x*=x)%=X,y>>=1;return t;}
I int gcd(CI x,CI y) {return y?gcd(y,x%y):x;}
I void Calc()//计算贡献
{
	RI i,j,S=Fac[n],T=0;for(i=1;i<=t;++i) (T+=c[i]*(l[i]/2))%=X-1;//T统计总长度：同一循环
	for(i=1;i<=t;++i) for((T+=c[i]*(c[i]-1)/2*l[i])%=X-1,j=i+1;j<=t;++j) (T+=c[i]*c[j]*g[l[i]][l[j]])%=X-1;//T统计总长度：不同循环
	for(i=1;i<=t;++i) (S*=QP(Inv[l[i]],c[i])*IFac[c[i]])%=X;(ans+=S*pw[T])%=X;//S统计方案数，ans统计答案
}
I void dfs(CI x,CI s)//搜索划分方案
{
	if(s==n) return Calc();if(x==1) return (void)(l[++t]=1,c[t]=n-s,Calc(),--t);//边界
	dfs(x-1,s),++t;for(RI i=1;x*i<=n-s;++i) l[t]=x,c[t]=i,dfs(x-1,s+x*i);--t;//枚举C
}
int main()
{
	RI i,j;for(scanf("%d",&n),pw[0]=i=1;i^X;++i) pw[i]=(pw[i-1]<<1)%X;//预处理2的幂
	for(Fac[0]=i=1;i<=n;++i) Fac[i]=Fac[i-1]*i%X;//预处理阶乘
	for(IFac[n]=QP(Fac[n],X-2),i=n-1;~i;--i) IFac[i]=IFac[i+1]*(i+1)%X;//预处理阶乘逆元
	for(i=1;i<=n;++i) for(Inv[i]=QP(i,X-2),j=1;j<=n;++j) g[i][j]=gcd(i,j);//预处理逆元和gcd
	return dfs(n,0),printf("%d\n",IFac[n]*ans%X),0;//dfs后输出答案
}