polya定理
算法简述
Burnside定理：设G是上的置换群，G是N上可引出不同的等价类，其中不同的等价类的个数为,其中，为置换g中不变元的个数，及g中1阶循环的个数。
polya定理：设是n个对象的置换群，用m种颜色给这n个对象着色，则不同的着色方案数为：。其中为置换的循环节数。
主要运用场合及思路
Burnside定理和Polya定理是组合数学中，用来计算全部互异的组合状态的个数的一个十分高效、简便的工具。一般在题目中出现经过旋转、翻转等方式重合的计数问题时往往会用到polya定理或Burnside定理，但由于题目中的取值范围很大，有时就会用到欧拉函数进行优化。
它们求解的一般步骤是：
- 确定置换群。
- polya：计算每个置换下循环节个数。Burnside：求解每个置换下本质不同的方案数即在该置换下保持不变的方案数。有时会用到组合数学或动态规划的方法进行计数。
- 代入公式得到答案。
模版
- poj2409(裸的polya)
#include<cstdio> typedef long long LL; LL Pow(LL a, LL b) { LL ans; for (ans = 1; b; b >>= 1) { if (b & 1) ans *= a; a *= a; } return ans; } int GCD(int x, int y) { return y ? GCD(y, x % y) : x; } int main() { int n, k, i; LL ans; while (scanf("%d%d", &k, &n), n || k) { if (n & 1) ans = Pow(k, n / 2 + 1) * n; else ans = Pow(k, n / 2 + 1) * (n / 2) + Pow(k, n / 2) * (n / 2); for (i = 1; i <= n; i++) ans += Pow(k, GCD(n, i)); ans = printf("%lld\n", ans / (2 * n)); } return 0; } - 用欧拉函数进行优化，可用dfs或循环的方式枚举因子，一般形式如下：
//循环枚举因子。 for(int i = 1;i*i <= n;i ++){ if(i * i == n){ ans = (ans + pow(n,i-1,p)*eular(i,p)%p) % p; } else if(n % i == 0){ ans += (pow(n,i-1,p)*eular(n/i,p)%p + pow(n,n/i-1,p)*eular(i,p)%p) % p; ans %= p; } } dfs枚举因子。cntf为素因子个数，num数组存放素因子 void dfs(int i, int now) { if(i == cntf) { ans = (ans + eular(now)*pow(n, n/now - 1))%MOD; if(now*now != n) ans = (ans + eular(n/now)*pow(n, now - 1))%MOD; return ; } int tmp = 1, j; for(j = 0; j <= num[i]; ++j) { if(LL(now*tmp)*LL(now*tmp) > n) return ; dfs(i + 1, now*tmp); tmp *= fac[i]; } }
例题
hdu5080 Colorful Toy
题意：平面上有n个点m条边，用c种颜色给这m条边染色，如果两个方案经过旋转后重合则视为一种方案，求有多少种染色方案。
思路：此题就是简单的polya定理的应用，关键在于枚举旋转角度后对应点的处理。
参考代码：
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <vector>
using namespace std;
#include <math.h>
#define LL long long
#define EXIT exit(0);
#define DEBUG puts("Here is a BUG");
#define CLEAR(name, init) memset(name, init, sizeof(name))
const double eps = 1e-8;
const int MAXN = (int)1e9 + 5;
using namespace std;
const LL mod = 1LL*1000000007;
#define Vector Point
int dcmp(double x)
{
return fabs(x) < eps ? 0 : (x < 0 ? -1 : 1);
}
struct Point
{
double x, y;
Point(const Point& rhs): x(rhs.x), y(rhs.y) { } //拷贝构造函数
Point(double x = 0.0, double y = 0.0): x(x), y(y) { } //构造函数
friend istream& operator >> (istream& in, Point& P)
{
return in >> P.x >> P.y;
}
friend ostream& operator << (ostream& out, const Point& P)
{
return out << P.x << ' ' << P.y;
}
friend Vector operator + (const Vector& A, const Vector& B)
{
return Vector(A.x+B.x, A.y+B.y);
}
friend Vector operator - (const Point& A, const Point& B)
{
return Vector(A.x-B.x, A.y-B.y);
}
friend Vector operator * (const Vector& A, const double& p)
{
return Vector(A.x*p, A.y*p);
}
friend Vector operator / (const Vector& A, const double& p)
{
return Vector(A.x/p, A.y/p);
}
friend bool operator == (const Point& A, const Point& B)
{
return dcmp(A.x-B.x) == 0 && dcmp(A.y-B.y) == 0;
}
friend bool operator < (const Point& A, const Point& B)
{
return A.x < B.x || (A.x == B.x && A.y < B.y);
}
friend bool operator != (const Point &A,const Point &B)
{
return dcmp(A.x-B.x) != 0 || dcmp(A.y-B.y) != 0;
}
void in(void)
{
scanf("%lf%lf", &x, &y);
}
void out(void)
{
printf("%lf %lf", x, y);
}
}p[55],p2[55];
double Dot(const Vector& A, const Vector& B)
{
return A.x*B.x + A.y*B.y; //点积
}
double Length(const Vector& A)
{
return sqrt(Dot(A, A)+eps);
}
double Angle(const Vector& A, const Vector& B)
{
return acos(Dot(A, B)/Length(A)/Length(B)); //向量夹角
}
double Cross(const Vector& A, const Vector& B)
{
return A.x*B.y - A.y*B.x; //叉积
}
//向量绕起点旋转
Vector Rotate(const Vector& A, const double& rad)
{
return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}
bool cmp(int A,int B)
{
if(dcmp(Cross(p[A],p[B]) != 0))
{
return Cross(p[A],p[B]) > 0;
}
else return Length(p[A]) < Length(p[B]);
}
int n,m,c;
int Next[55];
bool flag;
vector <int> vec;
int map[55][55];
bool check1(double ang,int dif){
int sz = vec.size();
for(int i = 0;i < sz;i ++){
p2[Next[vec[i]]] = Rotate(p[vec[i]],ang);
if(p2[Next[vec[i]]] != p[Next[vec[i]]]) return false;
}
return true;
}
bool check2(int dif){
for(int i = 0;i < n;i ++)
for(int j = 0;j < n;j ++)
if(map[i][j] != map[Next[i]][Next[j]]) return false;
return true;
}
long long quick_mod(LL a,int k){
LL ans = 1;
while(k){
if(k&1){
ans *= a;
ans %= mod;
}
a = a * a;
a %= mod;
k >>= 1;
}
return ans;
}
bool vis[100] ;
long long cal(int dif){
memset(vis,false,sizeof(vis));
int num = 0;
for(int i = 0;i < n;i ++){
if(!vis[i]){
vis[i] = true;
num ++;
int j = Next[i];
while(j != i){
vis[j] = true;
j = Next[j];
}
}
}
return quick_mod(1LL*c,num);
}
void gcd(LL a,LL b,LL& d,LL & x,LL & y){
if(!b) { d = a ; x = 1;y = 0;}
else {gcd(b,a%b,d,y,x); y-=x*(a/b);}
}
LL inv(LL a){
LL d,x,y;
gcd(a,mod,d,x,y);
return d == 1 ? (x + mod) % mod : -1;
}
int main(){
int t;
// freopen("a.txt","r",stdin);
scanf("%d",&t);
while(t --){
scanf("%d%d%d",&n,&m,&c);
double sumx = 0,sumy = 0;
vec.clear();
memset(map,0,sizeof(map));
for(int i = 0;i < n;i ++){
p[i].in();
sumx += p[i].x;
sumy += p[i].y;
}
Point center(sumx/n,sumy/n);
for(int i = 0;i < n;i ++){
p[i] = p[i]-center;
}
flag = false;
center.x = 0;
center.y = 0;
for(int j = 0;j < n;j ++){
if(p[j] == center) Next[j] = j;
else vec.push_back(j);
}
for(int i = 0;i < m;i ++){
int u,v;
scanf("%d%d",&u,&v);
u--;v--;
map[u][v] = map[v][u] = 1;
}
int ti = 0;
long long ans = 0;
int sz = vec.size();
sort(vec.begin(),vec.end(),cmp);
for(int i = 0;i < sz;i ++){
for(int j = 0;j < sz;j ++){
Next[vec[j]] = vec[(j+i)%sz];
}
Vector tmp(p[vec[0]]);
Vector tmp2(p[Next[vec[0]]]);
double ang = atan2(tmp2.y,tmp2.x) - atan2(tmp.y,tmp.x);
int dif = i;//
if(check1(ang,dif) && check2(dif)){
ans += cal(dif);
ans %= mod;
ti ++;
}
}
ans *= inv(ti);
ans %= mod;
printf("%lld\n",ans);
}
return 0;
}
bzoj1488
题意：求两两互不同构的含n个点的简单图有多少种。
思路：点重新标号的过程就是置换的过程，对于两点之间有无边相连可以想象成有每条边可以染无色和有色两种颜色，那么我们就需要求出每种置换对应边循环的个数。经过分析我们可以知道每一种共轭类对应的置换其边置换是一样的。然后我们考虑边置换的循环节个数，当一条边的两个端点在一个长度为的循环中时会产生个边循环节，当一条边的两个端点分别在两个长度为的循环中时会产生个边循环。对于形如的共轭类，会有个元素，我们可dfs枚举每一种共轭类便可求解。
参考代码：
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
#define N 65
#define mod 997
int fac[N];
int n;
int quick(int n,int k){
int ans = 1;
while(k){
if(k & 1) ans = ans * n % mod;
n = n * n % mod;
k >>= 1;
}
return ans;
}
int _inv(int n){
return quick(n,mod-2);
}
void ini(){
fac[0] = 1;
for(int i = 1;i < N;i ++){
fac[i] = fac[i-1] * i % mod;
}
}
int num[N][2];
int cnt,ans;
int gcd(int n,int m){
return m == 0 ? n : gcd(m,n % m);
}
void cal(){
int x = 0;
for(int i = 0;i < cnt;i ++){
x += num[i][1]*(num[i][1]-1)/2 * num[i][0] + num[i][0] / 2 * num[i][1];
for(int j = i+1;j < cnt;j ++){
x += num[i][1] * num[j][1] * gcd(num[i][0],num[j][0]);
}
}
int res = quick(2,x);
res = res * fac[n] % mod;
for(int i = 0;i < cnt;i ++){
res = res* _inv(fac[num[i][1]]) % mod * _inv(quick(num[i][0],num[i][1])) % mod;
res %= mod;
}
ans = (ans + res) % mod;
}
void dfs(int now,int left){
if(left == 0){
cal();
return ;
}
if(now > left) return ;
dfs(now + 1,left);
for(int i = 1;i * now <= left;i ++){
num[cnt][0] = now;
num[cnt++][1] = i;
dfs(now+1,left-i*now);
cnt --;
}
}
int main(){
ini();
while(~scanf("%d",&n)){
cnt = 0,ans = 0;
dfs(1,n);
ans = ans * _inv(fac[n]) % mod;
printf("%d\n",ans);
}
return 0;
}
poj2154 Color
题意：n种颜色对一串n个珠子的项链进行染色，求有多少种染色方案(若两条项链经过旋转后能完全重合，则视为同种方案)(n<=1000000000)。
思路：此题为经典的polya题目，将项链顺时针旋转i格后，其循环节个数为，所以此题的最终结果应为。但是n的数目为1000000000，枚举i会超时。我们不妨换个角度想，的个数是很少的，所以我们可以枚举每个然后乘上相应的个数就是答案了。假设，所以的个数也就是所有不大于的数中与其互质的数目，我们可以用欧拉函数或容斥原理求得。
参考代码：
#include <iostream>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
#define N 36000
int pri[N];
bool is[N];
int pcnt;
void ini(){
pcnt = 0;
memset(is,false,sizeof(is));
for(int i = 2;i < N;i ++){
if(!is[i]){
pri[pcnt ++] = i;
for(int j = i+i;j < N;j += i)
is[j] = true;
}
}
}
int eular(int n,int p){
int ans = n;
for(int i = 0;i < pcnt && pri[i] * pri[i] <= n;i ++){
if(n % pri[i] == 0){
ans = ans - ans / pri[i];
while(n % pri[i] == 0)
n /= pri[i];
}
}
if(n > 1) ans = ans - ans/n;
return ans % p;
}
int quick(int n,int k,int p){
int ans = 1;
n %= p;
while(k){
if(k&1) ans = ans * n % p;
n = n * n % p;
k >>= 1;
}
return ans;
}
int main(){
ini();
int t;
scanf("%d",&t);
while(t --){
int n,p;
scanf("%d%d",&n,&p);
int ans = 0;
for(int i = 1;i*i <= n;i ++){
if(i * i == n){
ans = (ans + quick(n,i-1,p)*eular(i,p)%p) % p;
}
else if(n % i == 0){
ans += (quick(n,i-1,p)*eular(n/i,p)%p + quick(n,n/i-1,p)*eular(i,p)%p) % p;
ans %= p;
}
}
printf("%d\n",ans);
}
return 0;
}
POJ 2888 Magic Bracelet
题意：给一条n个珠子的项链用m种颜色进行染色，其中有k对颜色不能出现在相邻的珠子上，如果两种方案可通过旋转使其一致则视为一种方案，求有多少种不同方案数。
思路：此题可用Burnside定理做，一共有n种置换，对于每种置换我们需要求出在该置换下保持不变的着色方案数。由于项链每旋转i个珠子，就会有个长度为的循环节，而且可以发现项链中任意相邻的长度为的珠子处于不同的循环中，所以此时保持不变的着色方案数等价于长度为的项链满足合法条件的方案数。构造一个合法矩阵g，其中表示颜色可以相邻，表示不可以，将其看作一个无向图，那么表示一个点经过k条路之后到达某点的方法数，此时表示从种颜色开始染色最后染完一条项链的方案数，将其累加就是最后的答案。
参考代码：
#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
const int max_n=15;
const int MOD=9973;
int N;
struct Mat
{
int mat[max_n][max_n];
};
Mat operator*(Mat a,Mat b)
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
for(int k=0; k<N; k++)
{
if(a.mat[i][k] && b.mat[k][j])
{
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%MOD;
}
}
}
}
return c;
}
Mat operator^(Mat A,int x)
{
Mat c;
memset(c.mat,0,sizeof(c.mat));
for(int i=0; i<N; i++) c.mat[i][i]=1;
for(; x; x>>=1)
{
if(x&1)
{
c=c*A;
}
A=A*A;
}
return c;
}
int mypow(int a,int b)
{
int res=1;
for(; b; b>>=1)
{
if(b&1)
{
res*=a;
res%=MOD;
}
a*=a;
a%=MOD;
}
return res;
}
int M,K;
int PHI(int n)
{
int i,res=n;
long long j;
for(i=2,j=4LL; j<=(long long)n; i++,j+=i+i-1)
{
if(!(n%i))
{
res=res/i*(i-1);
while(!(n%i)) n/=i;
}
}
if(n>1) res=res/n*(n-1);
return res%MOD;
}
Mat A;
int solve(int p)
{
int res=0;
Mat C=A^p;
for(int i=0; i<N; i++)
{
res+=C.mat[i][i];
res%=MOD;
}
return res;
}
int main()
{
int T;
scanf("%d",&T);
for(int ncas=1; ncas<=T; ncas++)
{
scanf("%d%d%d",&M,&N,&K);
int u,v;
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
A.mat[i][j]=1;
}
}
for(int i=0; i<K; i++)
{
scanf("%d%d",&u,&v);
A.mat[u-1][v-1]=A.mat[v-1][u-1]=0;
}
int ans=0;
for(int p=1; p*p<=M; p++)
{
if(M%p==0)
{
if(p*p==M)
{
ans = (ans+ PHI(p)*solve(p) )%MOD;
}
else
{
ans = (ans+ PHI(p)*solve(M/p)+ PHI(M/p)*solve(p) )%MOD;
}
}
}
int inv=mypow((M%MOD),MOD-2);
ans*=inv;
ans%=MOD;
printf("%d\n",ans);
}
return 0;
}