BZOJ 1488 Luogu P4727 [HNOI2009]图的同构 (Burnside引理、组合计数)

题目链接

(Luogu) https://www.luogu.org/problem/P4727
(BZOJ) https://www.lydsy.com/JudgeOnline/problem.php?id=1488

题解

Burnside引理经典题。

首先考虑一个\(O(n!\times poly(n))\)暴力: 枚举点的置换，然后计算在置换下保持不变的图的个数。
把置换拆成若干个轮换。
(1) 考虑轮换内部: 假设一轮换为\((a_1\ a_2\ ...\ a_n)\), 那么\((a_1,a_2),(a_2,a_3),...,(a_n,a_1)\)这些边要么都存在要么都不存在；\((a_1,a_3),(a_2,a_4),...,(a_{n-1},a_{1}),(a_{n},a_2)\)这些边也要么都存在要么都不存在；一般地说，对于任何一个\(d\), 所有的\((a_i,a_{(i+d)\mod n})\)这些边要么都存在要么都不存在，因此轮换内部一共有\(2^{\frac{n}{2}}\)种方案。
(2) 考虑轮换之间: 假设两轮换分别为\((a_1,a_2,...,a_n),(b_1,b_2,...,b_m)\)则有: \((a_1,b_1),(a_2,b_2),...(a_i,b_i)\)这些边存在情况都相同；\((a_1,b_2),(a_2,b_3),...,(a_i,b_{i+1})\)这些边存在情况都相同；以此类推，可以得到两轮换之间共有\(2^{\gcd(n,m)}\)种方案。
所有的置换方案数相加，最后除以置换总数\(n!\).

然后现在考虑\(n\le 60\)怎么办。
当\(n\le 60\)时，我们可以枚举拆分数(\(60\)的拆分数约为百万级别)。
已知一个拆分的方案(方案是指一个无标号序列\(a\)满足\(\sum a_i=n\)，其长度为\(cnt\))，它对应了多少个不同排列的轮换分拆？
首先，长度为\(L\)的轮换共\((L-1)!\)种。
然后我们要处理标号问题。
假设轮换之间是有区别的，那么标号方案数为\(\frac{n!}{\prod^{cnt}_{i=1}a_i!}\).
但是长度相等的轮换之间没有区别，所以除以\(\prod {num_i!}\), 其中\(num_i\)表示\(i\)在\(a\)中的出现次数。
最后乘起来得到\(\frac{n!}{\prod^{cnt}_{i=1}a_i\prod^n_{i=1}num_i!}\)
累加即可。

时间复杂度？
枚举所有拆分方案，求\(cnt^2\)之和，我用程序计算得当\(n=60\)时该值约为\(2.7\times 10^8\).

代码

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cassert>
#include<iostream>
#define llong long long
using namespace std;

inline int read()
{
	int x=0; bool f=1; char c=getchar();
	for(;!isdigit(c);c=getchar()) if(c=='-') f=0;
	for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0');
	if(f) return x;
	return -x;
}

const int N = 60;
const int P = 997;
llong fact[N+3],finv[N+3],inv[N+3];
llong pw2[N+3];
int a[N+3];
int num[N+3];
int gcd[N+3][N+3];
int n,cnt;
llong ans;

int GCD(int x,int y)
{
	return y==0?x:GCD(y,x%y);
}

llong quickpow(llong x,llong y)
{
	llong cur = x,ret = 1ll;
	for(int i=0; y; i++)
	{
		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}
		cur = cur*cur%P;
	}
	return ret;
}

llong calc()
{
	llong ret = fact[n];
	for(int i=1; i<=cnt; i++)
	{
		ret = ret*inv[a[i]]%P;
	}
	for(int i=1; i<=n; i++)
	{
		ret = ret*finv[num[i]]%P;
	}
	for(int i=1; i<=cnt; i++)
	{
		ret = ret*pw2[a[i]>>1]%P;
	}
	for(int i=1; i<=cnt; i++)
	{
		for(int j=i+1; j<=cnt; j++)
		{
			ret = ret*pw2[gcd[a[i]][a[j]]]%P;
		}
	}
	return ret;
}

void dfs(int sum)
{
	if(sum==n)
	{
		ans = (ans+calc())%P;
		return;
	}
	for(int i=a[cnt]; i+sum<=n; i++)
	{
		cnt++; a[cnt] = i; num[i]++;
		dfs(i+sum);
		a[cnt] = 0; cnt--; num[i]--;
	}
}

int main()
{
	pw2[0] = 1ll; for(int i=1; i<=N; i++) pw2[i] = (pw2[i-1]<<1)%P;
	fact[0] = 1ll; for(int i=1; i<=N; i++) fact[i] = fact[i-1]*i%P;
	finv[N] = quickpow(fact[N],P-2); for(int i=N-1; i>=0; i--) finv[i] = finv[i+1]*(i+1)%P;
	for(int i=1; i<=N; i++) inv[i] = finv[i]*fact[i-1]%P;
	for(int i=1; i<=N; i++) for(int j=1; j<=N; j++) gcd[i][j] = GCD(i,j);
	scanf("%d",&n);
	if(n==0) {printf("1"); return 0;}
	a[0] = 1; dfs(0);
	ans = ans*finv[n]%P;
	printf("%lld\n",ans);
	return 0;
}