[POJ 2888]Magic Bracelet[Polya(Burnside) 置换 矩阵]

也许刚好的阅读体验
$Description$Description

大意：给一条长度为 $n$n的项链，有 $m$m种颜色，另有 $k$k条限制，每条限制为不允许 $x,y$x,y颜色连在一起。要求有多少种本质不同的染色方式，本质不同的两种染色方式必须旋转不能互相得到。
输入方式：
第一行 $t,$t,表示t组数据
接下来 $t$t组数据：
每组数据第一行为 $n,m,k$n,m,k
接下来 $k$k行，每行两个数 $x,y$x,y表示不允许 $x,y$x,y颜色连在一起。
答案对9973取模
$(1\le n\le 10^9,gcd(n,9973)=1),m(1\le m\le 10),k(1\le k\le m(m-1\left)/2\right)$(1≤n≤109,gcd(n,9973)=1),m(1≤m≤10),k(1≤k≤m(m−1)/2)
$Solution$Solution

本篇题解假设大家都会没有 $k$k条限制的版本

由于染色有限制，所以 $polya$polya定理就不好用了
用 $burnside$burnside定理解决(发现大家标题都打得polya,于是也这么打标题了)

首先枚举不同的置换，即枚举循环长度，这一块和用 $polya$polya定理有点像
由于 $n$n很大，要用 $\phi$φ来优化
一个置换中，循环的元素的所有颜色必须相同，所以我们要计算有多少个循环，这些循环有多少种染色方法

计算染色方法
考虑在该置换内一个循环一个循环的染色，我们可以看做从不同颜色节点走向另一节点
用邻接矩阵 $f\left[i\right]\left[j\right]$f[i][j]表示从 $i$i能否走向 $j$j
除去那 $k$k条限制，所有其他的 $i,j$i,j , $f\left[i\right]\left[j\right]=1$f[i][j]=1,即在染为第 $i$i种颜色后是可以染为第 $j$j种颜色的

这时候想到邻接矩阵的妙用，用矩乘计算 $n$n次之后得到的 $f\left[i\right]\left[j\right]$f[i][j]的意义发生改变
$f\left[i\right]\left[j\right]$f[i][j]表示从 $i$i出发，走 $n$n步，最后结束在 $j$j有多少种走法。
由于是一个走一圈，所以 从 $i$i出发应在 $i$i结束，答案贡献为 $f\left[i\right]\left[i\right]$f[i][i]
这样就可以计算长度为 $n$n的循环的染色方法了：求出原矩阵的 $n$n次方后的矩阵后 ${\sum }_{i=1}^{m}f\left[i\right]\left[i\right]$∑i=1m​f[i][i]

注意取模
代码

/******************************* Author:Morning_Glory LANG:C++ Created Time:2019年07月04日 星期四 14时25分13秒 *******************************/
#include <cstdio>
#include <fstream>
#include <cstring>
using namespace std;
const int maxn = 25;
const int mod = 9973;
//{{{cin 读入优化
struct IO{
	template<typename T>
	IO & operator>>(T&res){
		res=0;
		bool flag=false;
		char ch;
		while((ch=getchar())>'9'||ch<'0')	 flag|=ch=='-';
		while(ch>='0'&&ch<='9') res=(res<<1)+(res<<3)+(ch^'0'),ch=getchar();
		if (flag)	 res=~res+1;
		return *this;
	}
}cin;
//}}}
int t,n,m,k,x,y,ni,ans;
//{{{Matrix
struct Matrix{
	int rec[maxn][maxn];
	Matrix(){	memset(rec,0,sizeof(rec));}
	Matrix operator = (const Matrix &y){
		for (int i=1;i<=m;++i)
			for (int j=1;j<=m;++j)
				rec[i][j]=y.rec[i][j];
		return *this;
	}
	friend Matrix operator * (const Matrix &a,const Matrix &b){
		Matrix t;
		for (int i=1;i<=m;++i)
			for (int j=1;j<=m;++j)
				for (int k=1;k<=m;++k)
					t.rec[i][j]=(t.rec[i][j]+a.rec[i][k]*b.rec[k][j])%mod;
		return t;
	}
	Matrix operator ^ (int b){
		Matrix s,a;
		a=*this;
		for (int i=1;i<=m;++i)	s.rec[i][i]=1;
		for (;b;b>>=1,a=a*a)
			if (b&1)	s=s*a;
		return s;
	}
}a,b;
//}}}Martix
//{{{get_phi
int get_phi (int x)
{
	int res=x;
	for (int i=2;i*i<=x;++i)
		if (x%i==0){
			res=res/i*(i-1);
			while (x%i==0)	x/=i;
		}
	if (x>1)	res=res/x*(x-1);
	return res%mod;//此处要取模
}
//}}}
//{{{ksm
int ksm (int a,int b)
{
	int s=1;
	a%=mod;
	for (;b;a=1ll*a*a%mod,b>>=1)
		if (b&1)	s=1ll*s*a%mod;
	return s;
}
//}}}
//{{{calc
int calc (int x)
{
	b=a^x;
	int res=0;
	for (int i=1;i<=m;++i)	res=(res+b.rec[i][i])%mod;
	return res;
}
//}}}
int main()
{
	cin>>t;
	while (t--){
		cin>>n>>m>>k;
		ans=0,ni=ksm(n,mod-2);//ni n的逆元
		for (int i=1;i<=m;++i)
			for (int j=1;j<=m;++j)
				a.rec[i][j]=1;
		while (k--){
			cin>>x>>y;
			a.rec[x][y]=a.rec[y][x]=0;
		}
		for (int i=1;i*i<=n;++i)
			if (n%i==0){
				ans=(ans+get_phi(i)*calc(n/i)%mod)%mod;
				if (i*i!=n)	ans=(ans+get_phi(n/i)*calc(i)%mod)%mod;
			}
		ans=ans*ni%mod;
		printf("%d\n",ans);
	}
	return 0;
}