首页 > 代码库 > HDU 3923 Invoker 【裸Polya 定理】
HDU 3923 Invoker 【裸Polya 定理】
参考了 http://blog.csdn.net/ACM_cxlove?viewmode=contents by---cxlove 的模板
对于每一种染色，都有一个等价群，例如旋转，翻转等。我们将每一种变换转换成一个置换群，通过置换群得到的都是等价的染色方案
最终我们要求的是非等价的染色方案数。
在Burnside定理中给出，在每一种置换群也就是等价群中的数量和除以置换群的数量，即非等价的着色数等于在置换群中的置换作用下保持不变的着色平均数。
我们以POJ 2409 Let it Bead为例http://poj.org/problem?id=2409
N个物品的环，M种颜色。非常基础的问题
那么总共有多少个置换群呢，显然旋转是等价的，那么包括旋转0度，总共便是N个旋转。
对于每一个旋转的等价数量怎么算呢。我们发现每一个旋转置换中有若干个循环节，例如1，3，5是一个循环节，那么旋转一次，1号位到了3号位，3号位到了5号位，显然一个循环节内的颜色一定是一样的。有M种颜色选择，问题转换成有多少个循环节。对于旋转显然有GCD(i,n)。那只要枚举i便可以，这就是Polya定理。
而在翻转当中，如果N为奇数，明显有一个物品是不动的，其它的两两对称，颜色也是一样的。n/2+1个循环节，
如果 N为偶数，分为两种情况，一种是对称轴不过物品，那么所有物品两两对应，n/2个循环节，另外一种是对称轴经过两个物品，(n-2)/2+1+1个循环节。对于每一个循环节有M种颜色可以选择。
#include<iostream>#include<cstring>#include<queue>#include<cstdio>#include<cmath>#include<algorithm>#define N 100005#define inf 1<<29#define MOD 2007#define LL long longusing namespace std;LL gcd(LL a,LL b){ return b==0?a:gcd(b,a%b);}LL Pow(LL a,int b){ LL ret=1; while(b){ if(b&1) ret=ret*a; a=a*a; b>>=1; } return ret;}LL Polya(int n,int m){ LL sum=0; //枚举n种旋转 for(int i=1;i<=n;i++) //每个循环节是m种颜色可选 //总共有gcd(n,i)个循环节 sum+=Pow(m,gcd(n,i)); if(n&1) //如果为奇数，所有位置上的循环节数量都为n/2+1 sum+=n*Pow(m,n/2+1); else //否则要分奇偶，各一半 sum+=n/2*Pow(m,n/2)+n/2*Pow(m,n/2+1); return sum/2/n;}int n,m;int main(){ while(scanf("%d%d",&m,&n)!=EOF&&n+m){ printf("%lld\n",Polya(n,m)); } return 0;}
对于N如果非常大的话，枚举旋转就非常耗时。接下来可以有个优化
我们枚举循环节长度L，那么循环节个数便是N/L,必然L必须是N的约数才行。
那么d=N/L=gcd(i,N)，对于每一个L，我们需要求出有多少个i满足左边的式子。
令i=d*t,gcd(d*t,L*d)=d，要左边的式子成立，明显gcd(t,L)==1,否则最大公约数不为d。
那么对于任意的t满足与L互质即可，便是L的欧拉函数值
这样就可以在sqrt(n)的复杂度内枚举L， ∑(phi(L) * M^(N/L) ) % p (L即枚举值) 。
#include<iostream>#include<cstring>#include<queue>#include<cstdio>#include<cmath>#include<algorithm>#define N 100005#define inf 1<<29#define MOD 2007#define LL long longusing namespace std;LL gcd(LL a,LL b){ return b==0?a:gcd(b,a%b);}LL Eular(LL n){ LL ret=1; for(int i=2;i*i<=n;i++){ if(n%i==0){ n/=i; ret*=i-1; while(n%i==0){ n/=i; ret*=i; } } } if(n>1) ret*=n-1; return ret;}LL Pow(LL a,LL b){ LL ret=1; while(b){ if(b&1) ret=ret*a; a=a*a; b>>=1; } return ret;}LL Polya(int n,int m){ LL sum=0; int i; for(i=1;i*i<n;i++) if(n%i==0){ sum+=Eular(i)*Pow(m,n/i); sum+=Eular(n/i)*Pow(m,i); } if(i*i==n) sum+=Eular(i)*Pow(m,i); if(n&1) sum+=n*Pow(m,n/2+1); else sum+=n/2*Pow(m,n/2)+n/2*Pow(m,n/2+1); return sum/2/n;}int n,m;int main(){ while(scanf("%d%d",&m,&n)!=EOF&&n+m) printf("%lld\n",Polya(n,m)); return 0;}
Source Code:
#include<iostream>#include<cstring>#include<cstdio>#include<cmath>#define MOD 1000000007 #define LL long longusing namespace std;int prime[30]={2,3,5,7,11,13,17,19,23,29,31,37,41,43, 47,53,59,61,67,71,73,79,83,89,97,101},cnt=25;LL n,m;LL eular(LL n){ LL sum = 1; for(int i = 2; i <= sqrt(1.0 + n); ++i) if(n % i==0){ sum *= (i-1); n /= i; while(n % i == 0){ sum *= i; n /= i; } } if(n > 1) sum *= (n-1); return sum % MOD;}LL Pow(LL a,LL b){ LL ret=1; while(b){ if(b&1) ret=(ret*a)%MOD; a=(a*a)%MOD; b>>=1; } return ret;}LL Polya(){ LL sum=0,i; for(i=1;i*i<n;i++){ if(n%i==0){ sum=(sum+eular(i)*Pow(m,n/i))%MOD; sum=(sum+eular(n/i)*Pow(m,i))%MOD; } } if(i*i==n) sum=(sum+eular(i)*Pow(m,i))%MOD; if(n&1) sum=(sum+n*Pow(m,n/2+1))%MOD; else sum=(sum+n/2*(Pow(m,n/2)+Pow(m,n/2+1)))%MOD; return (sum*Pow(2*n,MOD-2))%MOD;}int main(){ LL t,cas=0; scanf("%I64d",&t); while(t--){ scanf("%I64d%I64d",&m,&n); printf("Case #%I64d: %I64d\n",++cas,Polya()); } return 0;}
HDU 3923 Invoker 【裸Polya 定理】