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! imes 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!).

然后现在考虑(nle 60)怎么办。
当(nle 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_iprod^n_{i=1}num_i!})
累加即可。

时间复杂度？
枚举所有拆分方案，求(cnt^2)之和，我用程序计算得当(n=60)时该值约为(2.7 imes 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
",ans);
	return 0;
}