等价类计数常见的环形模型

Burnside引理

把等价变换的所有置换列出，构成置换群。此时特别注意群的封闭性、且不重复不遗漏地覆盖等价类的所有的置换，这是最难的部分。然后求出每个置换分解为轮换的个数(轨道的个数)，可以手推或用程序暴力计算(dfs或直接并查集合并连通分量)。那么本质不同的方案是群中每个置换下不变集合的 (染色时即 c^(每个置换的轨道数)) 的平均数。

等价类计数

变换下的等价类个数的计数问题，考虑使用Burnside引理消除等价类条件，以一定复杂度的代价转换为普通计数问题。

其中非常常见的是n元环旋转(有时还有翻转)。

仅旋转时，群的阶为n，易观察知其中第i个置换(旋转i个单位)的轨道数为gcd(i,n)。即置换的轨道数d均为n的因子(也是稳定子轨道定理的结论)。根据欧拉函数的定义，还可得出轨道数为d的置换个数为φ(n/d)个。

因此，旋转的问题中，可在O(sqrt(n))的时间内枚举n的因子，再对因子算欧拉函数得到置换的个数。这一部分直接暴力计算一般已经足够快，此外还可以继续优化。

若还可对环进行翻转，增加以n条对称轴为对称中心的置换即可。这些翻转的置换都非常简单，只需要注意奇偶不同以及经过轴上的点不动。

poj 1286 / poj 2409

基本的环旋转、翻转，裸的Polya定理。下为poj 1286代码。

```
/*
```
```
Author: fffasttime
```
```
Date:
```
```
 */
```
```
#include <iostream>
```
```
#include <cstdio>
```
```
#include <cmath>
```
```
#include <algorithm>
```
```
using namespace std;
```
```
typedef long long ll;
```
```
#define inc(i,n) for (int i=0;i<n;i++)
```
```
const int N=100010;
```
```
 
```
```
typedef long long ll;
```
```
 
```
```
ll po(ll a, ll x){
```
```
	ll ans=1;
```
```
	while (x--) ans*=a;
```
```
	return ans;
```
```
}
```
```
 
```
```
int main(){
```
```
	int n; while (cin>>n && n!=-1){
```
```
		if (n==0) {cout<<0<<'\n';continue;}
```
```
		ll sum=0;
```

		for (int i=0;i<n;i++) sum+=po(3,__gcd(i,n));

```
		if (!(n&1)){
```
```
			sum+=n/2*po(3,n/2+1);
```
```
			sum+=n/2*po(3,n/2);
```
```
		}
```
```
		else{
```
```
			sum+=n*po(3,n/2+1);
```
```
		}
```
```
		cout<<sum/n/2<<'\n';
```
```
	}
```
```
	return 0;
```
```
}
```

poj 2154 / luogu 【模板】Polya定理

就是环旋转的情况，使用上述的一点数论优化一下。

```
/*
```
```
Author: fffasttime
```
```
Date:
```
```
 */
```
```
#include <iostream>
```
```
#include <cstdio>
```
```
#include <cmath>
```
```
#include <algorithm>
```
```
using namespace std;
```
```
typedef long long ll;
```
```
#define inc(i,n) for (int i=0;i<n;i++)
```
```
const int N=100010;
```
```
 
```
```
ll a[110];
```
```
 
```
```
bool np[40000];
```
```
int pp[10000],ppc;
```
```
 
```
```
void init(){
```
```
	for (int i=2;i<40000;i++){
```
```
		if (!np[i]) pp[ppc++]=i;
```
```
		for (int j=i+i;j<40000;j+=i)
```
```
			np[j]=1;
```
```
	}
```
```
}
```
```
 
```
```
ll qpow(ll a, ll x, ll p){
```
```
	ll ans=1;
```
```
	for (;x;a=a*a%p,x>>=1)
```
```
		if (x&1) ans=ans*a%p;
```
```
	return ans;
```
```
}
```
```
 
```
```
int p[20],pc;
```
```
ll getphi(int x){
```
```
	ll ans=1;
```
```
	for (int i=0;i<pc;i++){
```
```
		if (p[i]*p[i]>x) break;
```
```
		if (x%p[i]==0){
```
```
			x/=p[i]; ans=ans*(p[i]-1);
```

			while (x%p[i]==0) x/=p[i],ans=ans*p[i];

```
		}
```
```
	}
```
```
	if (x>1) ans=ans*(x-1);
```
```
	return ans;
```
```
}
```
```
ll P;
```
```
int main(){
```
```
	int T; cin>>T;
```
```
	init();
```
```
	while (T--){
```

		int n0; cin>>n0>>P; int n=n0, r=sqrt(n+1);

```
		pc=0;
```
```
		for (int j=0;pp[j]<=r;j++){
```
```
			int i=pp[j];
```
```
			if (n%i==0) {
```
```
				p[pc++]=i;
```
```
				while (n%i==0) n/=i;
```
```
			}
```
```
		}
```
```
		if (n>1) p[pc++]=n;
```
```
		n=n0; ll ans=0;
```
```
		for (int i=1;i<=r;i++) if (n%i==0){
```
```
			ans+=getphi(n/i)*qpow(n,i-1,P);
```

			if (i*i<n) ans+=getphi(i)*qpow(n,n/i-1,P);

```
			ans%=P;
```
```
		}
```
```
		//cout<<ans<<'\n';
```
```
		cout<<ans<<'\n';
```
```
	}	
```
```
	return 0;
```
```
}
```

hdu 1812

也是Polya模板题，这题不是环旋转而是方阵旋转、翻转，不过也比较简单。旋转、翻转共8个置换，手推一下即可。

```
Huge qpow(Huge a, int x){
```
```
	Huge ans; ans.set("1");
```
```
	for (;x;a*=a,x>>=1)
```
```
		if (x&1) ans*=a;
```
```
	return ans;
```
```
}
```
```
 
```
```
int main(){
```
```
	int n;
```
```
	char c0[10];
```
```
	while (cin>>n>>c0){
```
```
		Huge sum; sum.set("0");
```
```
		Huge c; c.set(c0);
```
```
		Huge s2,s3; s2.set("2"); s3.set("3");
```
```
		sum+=qpow(c,n*n);
```
```
		if (n&1){
```
```
			sum+=s2*qpow(c,n*n/4+1);
```
```
			sum+=qpow(c,n*n/2+1);
```
```
			sum+=s2*s2*qpow(c,(n*n+n)/2);
```
```
		}
```
```
		else{
```
```
			sum+=s2*qpow(c,n*n/4);
```
```
			sum+=s3*qpow(c,n*n/2);
```
```
			sum+=s2*qpow(c,(n*n+n)/2);
```
```
		}
```
```
		sum.d2(); sum.d2(); sum.d2();
```
```
		cout<<sum<<'\n';
```
```
	}
```
```
	return 0;
```
```
}
```

poj 2888 带限制的等价类计数

同样是旋转变换视为等价。但这里有相邻不能取某些颜色的限制。因此不能再用Polya定理，只能直接用Burnside引理。根据Burnside引理，我们需要直接计算各个置换下不动的染色方案数。此时第i个旋转的轨道数仍为gcd(i,n)，轨道长度d=n/gcd(i,n)。考虑原先使用Polya定理时，c种染色不带限制，即长度为d的段可以任意染色，共c^d种方案。而此时可把长度为d的段取出，而段的首尾也需要满足限制条件，可视为长度为d的环。

于是问题转化为：对长度为d的环染c种颜色，某些颜色不能相邻，求方案数，旋转视为不同方案。这是一类比较经典的矩阵加速递推的计数问题，与https://www.luogu.org/problem/P1357类似。该递推还可使用路径计数解释：以颜色为节点，相邻合法的颜色连一条边。则邻接矩阵的k次幂即在图上走k步，选出长度为k的合法染色。长度为d的环即走d步后回到自身，即邻接矩阵d次方对角线的和。

这里枚举置换还是跟上面一样的套路，所以总复杂度乘上子问题的复杂度，O(sqrt(n) m^3 log(n))。

```
/*
```
```
Author: fffasttime
```
```
Date:
```
```
 */
```
```
#include <iostream>
```
```
#include <cstdio>
```
```
#include <cmath>
```
```
#include <algorithm>
```
```
#include <cstring>
```
```
using namespace std;
```
```
typedef long long ll;
```
```
#define inc(i,n) for (int i=0;i<n;i++)
```
```
const int N=100010;
```
```
 
```
```
const int P=9973;
```
```
 
```
```
int qpow(int a, int x){
```
```
	int ans=1;
```
```
	a%=P;
```
```
	for(;x;a=a*a%P,x>>=1)
```
```
		if (x&1) ans=ans*a%P;
```
```
	return ans;
```
```
}
```
```
 
```
```
int n,m,k;
```
```
 
```
```
struct Mat{
```
```
	int m[11][11];
```
```
	void cl(){memset(m,0,sizeof m);}
```
```
	void I(){cl();inc(i,n)m[i][i]=1;}
```
```
}b;
```
```
 
```

Mat operator*(const Mat &a, const Mat &b){

```
	Mat c; c.cl();
```
```
	inc(i,n) inc(j,n) inc(k,n)
```

		c.m[i][j]+=a.m[i][k]*b.m[k][j],c.m[i][j]%=P;

```
	return c;
```
```
}
```
```
Mat qpow(Mat a, int x){
```
```
	Mat ans; ans.I();
```
```
	for (;x;a=a*a,x>>=1)
```
```
		if (x&1) ans=ans*a;
```
```
	return ans;
```
```
}
```
```
int calc(int d){
```
```
	Mat a=b;
```
```
	Mat ans=qpow(a,d);
```
```
	int ret=0;
```
```
	inc(i,n) ret+=ans.m[i][i],ret%=P;
```
```
	return ret;
```
```
}
```
```
 
```
```
int phi(int n){
```
```
	int ans=n;
```
```
	for (int i=2;i*i<=n;i++)
```
```
		if (n%i==0){
```
```
			ans=ans/i*(i-1);
```
```
			while (n%i==0) n/=i;
```
```
		}
```
```
	if (n>1) ans=ans/n*(n-1);
```
```
	return ans%P;
```
```
}
```
```
 
```
```
int main(){
```
```
	int T; cin>>T;
```
```
	while (T--){
```
```
		cin>>m>>n>>k;
```
```
		inc(i,n) inc(j,n) b.m[i][j]=1;
```
```
		inc(i,k){
```
```
			int u,v; scanf("%d%d",&u,&v);
```
```
			u--; v--;
```
```
			b.m[u][v]=0;
```
```
			b.m[v][u]=0;
```
```
		}
```
```
		int ans=0;
```
```
		for (int i=1;i*i<=m;i++) if(m%i==0){
```
```
			ans+=phi(m/i)*calc(i);
```
```
			if (i*i!=m) ans+=phi(i)*calc(m/i);
```
```
			ans%=P;
```
```
		}
```
```
		printf("%d\n",ans*qpow(m,P-2)%P);
```
```
	}
```
```
	return 0;
```
```
}
```

hdu 2481 轮图生成树

这题除了增加了旋转等价的条件外，几乎和bzoj 1002轮状病毒一样，原问题递推为a[i]=3*a[i-1]-a[i-2]+2。有了上一题的套路，把bzoj 1002作为子问题即可。注意这里没有逆元，必须使用a/b%c=a%bc/b的方式计算。

```
/*
```
```
Author: fffasttime
```
```
Date:
```
```
 */
```
```
#include <iostream>
```
```
#include <cstdio>
```
```
#include <cmath>
```
```
#include <algorithm>
```
```
#include <cstring>
```
```
using namespace std;
```
```
typedef long long ll;
```
```
#define inc(i,n) for (int i=0;i<n;i++)
```
```
const int N=100010;
```
```
 
```
```
ll P=100000;
```
```
ll qmul(ll x,ll y){
```

	ll t=(x*y-(ll)((long double)x/P*y+0.5)*P);

```
	return t<0 ? t+P : t;
```
```
}
```
```
 
```
```
ll qpow(ll a, ll x){
```
```
	ll ans=1;
```
```
	a%=P;
```
```
	for(;x;a=qmul(a,a),x>>=1)
```
```
		if (x&1) ans=qmul(ans,a);
```
```
	return ans;
```
```
}
```
```
 
```
```
ll n;
```
```
 
```
```
struct Mat{
```
```
	ll m[3][3];
```
```
	void cl(){memset(m,0,sizeof m);}
```
```
	void I(){cl();inc(i,3)m[i][i]=1;}
```
```
};
```
```
Mat b={{{3,-1,2},{1,0,0},{0,0,1}}};
```
```
 
```

Mat operator*(const Mat &a, const Mat &b){

```
	Mat c; c.cl();
```
```
	inc(i,3) inc(j,3) inc(k,3)
```

		c.m[i][j]+=qmul(a.m[i][k],b.m[k][j]),c.m[i][j]%=P;

```
	return c;
```
```
}
```
```
Mat qpow(Mat a, ll x){
```
```
	Mat ans; ans.I();
```
```
	for (;x;a=a*a,x>>=1)
```
```
		if (x&1) ans=ans*a;
```
```
	return ans;
```
```
}
```
```
ll calc(int d){
```
```
	d--;
```
```
	Mat ans=qpow(b,d);
```
```
	ll ret=ans.m[0][0]+ans.m[0][2];
```
```
	return ret%P;
```
```
}
```
```
ll phi(int n){
```
```
	ll ans=n;
```
```
	for (ll i=2;i*i<=n;i++)
```
```
		if (n%i==0){
```
```
			ans=ans/i*(i-1);
```
```
			while (n%i==0) n/=i;
```
```
		}
```
```
	if (n>1) ans=ans/n*(n-1);
```
```
	return ans%P;
```
```
}
```
```
 
```
```
int main(){
```

	//for (int i=1;i<10;i++) cout<<calc(i)<<'\n';

```
	while (cin>>n>>P){
```
```
		P*=n;
```
```
		ll ans=0;
```
```
		for (int i=1;i*i<=n;i++) if(n%i==0){
```
```
			ans+=qmul(phi(n/i),calc(i));
```

			if (i*i!=n) ans+=qmul(phi(i),calc(n/i));

```
			ans%=P;
```
```
		}
```
```
		//cout<<ans<<'\n';
```
```
		printf("%lld\n",ans/n);
```
```
	}
```
```
	return 0;
```
```
}
```

spoj 422 矩阵转置交换次数

这题最大的困难之处在于它的置换特征非常不明显，需要仔细分析。由于矩阵长宽为2的幂次，首先寻找转置前后数组下标的二进制表示有什么特点。实际上是把当前下标二进制表示下a+b前后两段互换。而此交换可视为把一个01环向左循环旋转a个单位。我们需要求的是最少交换次数，对矩阵而言，一个01串可视为对矩阵一个元素的变换。由于一个轮换r需要r-1次交换操作才能完成，我们要求的答案就是总01串个数-轮换个数(即有多少个循环)。

然后注意到这样的轮换(轨道)个数即在左旋转a意义下等价的01串个数，也就是Polya定理要求的东西。想到这一点后就和上面的模板题差不太多了。这里旋转群的旋转间隔为gcd(a,a+b)，旋转置换群的阶为(a+b)/gcd(a,a+b)。唯一卡人的地方是这里n只有1e6而询问非常多，需要nlogn筛出n所有的因子，O(d(n))回答，d为因子个数函数。1e6内的d(n)不超过240，而平均只有10~20，因此不会超时。

```
#include <bits/stdc++.h>
```
```
using namespace std;
```
```
typedef long long ll;
```
```
 
```
```
const int P=1000003;
```
```
 
```
```
ll qpow(ll a, ll x){
```
```
	ll ans=1;
```
```
	for (;x;a=a*a%P,x>>=1)
```
```
		if (x&1) ans=ans*a%P;
```
```
	return ans;
```
```
}
```
```
 
```
```
const int maxn=1000010;
```
```
bool isp[maxn];
```
```
int p[maxn],phi[maxn],pc,pow2[maxn];
```
```
int fac[15000010],nxt[15000010],facc;
```
```
int head[maxn];
```
```
void add(int x, int d){
```
```
	nxt[facc]=head[x];
```
```
	fac[facc]=d;
```
```
	head[x]=facc++;
```
```
}
```
```
 
```
```
void gennumfunc(){
```
```
	memset(isp,1,sizeof(isp));
```
```
	isp[1]=0;  phi[1]=1;
```
```
	for (int i=2;i<maxn;i++){
```
```
		if (isp[i]) p[pc++]=i,phi[i]=i-1;
```

		for (int j=0;j<pc && i*p[j]<maxn;j++){

```
			isp[i*p[j]]=0;
```
```
			if (i%p[j]==0){
```
```
				phi[i*p[j]]=phi[i]*p[j];
```
```
				break;
```
```
			}
```

			phi[i*p[j]]=phi[i]*(p[j]-1); //f(ab)=f(a)f(b), when (a,b)=1

```
		}
```
```
	}
```
```
	pow2[0]=1;
```

	for (int i=1;i<maxn;i++) pow2[i]=pow2[i-1]*2%P;

```
	for (int i=1;i<maxn;i++)
```
```
		for (int j=i;j<maxn;j+=i)
```
```
			add(j,i);
```
```
}
```
```
int main(){
```
```
	gennumfunc();
```
```
	//cout<<clock()<<'\n';
```
```
	int T; cin>>T;
```
```
	while (T--){
```
```
		ll a,b; scanf("%lld%lld",&a,&b); 
```
```
		if (a*b==0){
```
```
			puts("0");
```
```
			continue;
```
```
		}
```
```
		ll g=__gcd(a,a+b);
```
```
		ll n=(a+b)/g;
```
```
		ll ans=0;
```
```
		for (int e=head[n];e;e=nxt[e]){
```
```
			int d=fac[e];
```
```
			ans+=1ll*pow2[d*g]*phi[n/d]%P;
```
```
			ans%=P;
```
```
		}
```

		printf("%lld\n",(qpow(2,n*g)+P-ans*qpow(n,P-2)%P)%P);

```
	}	
```
```
	return 0;
```
```
}
```