poj 2888 Magic Bracelet <polya定理>

题目：http://poj.org/problem?id=2888

题意：给定n（n <= 10^9）颗珠子，组成一串项链，每颗珠子可以用m种颜色中一种来涂色，如果两种涂色方法通过旋转项链可以得到视为等价。

　　　然后再给定K组限制，每组限制a、b代表颜色a和颜色b不能涂在相邻的珠子上面。问一共有多少种涂色方法。

思路：如果这题没有后面的限制，就和 poj 2154 一样了：http://www.cnblogs.com/jian1573/p/3234627.html

　　现在我们要处理的就是 K 种限制，可以用DP求解。 i为珠子编号， c为颜色编号那么：dp[i][c]=∑dp[i-1][cc] cc 为可以与 c 相邻的颜色编号；

　　由于N为 1e9 O(N) 会TLE，所以我们可以用矩阵快速幂来优化为O(lgN)：具体用m[i][j]=1,表示合法， m[i][j]=0表示不合法，

　　那么m^k后的对角线上的元素和即为所求。

  1 #include <iostream>
  2 #include <cmath>
  3 #include <cstdio>
  4 #include <cstring>
  5 using namespace std;
  6 const int Mod=9973;
  7 int N, M, K, T;
  8 struct Mar
  9 {
 10     int m[15][15];
 11     inline void zero( ){
 12         memset(m, 0, sizeof m);
 13     }
 14     inline void one( ){
 15         zero( );
 16         for( int i=0; i<M; ++i ){
 17             for( int j=0; j<M; ++j ){
 18                 m[i][j]=1;
 19             }
 20         }
 21     }
 22     inline void unit( ){
 23         zero();
 24         for( int i=0; i<M; ++ i )
 25             m[i][i]=1;
 26     }
 27     inline Mar operator *(const Mar &a) const {
 28         Mar C;C.zero();
 29         for( int i=0; i<M; ++i ){
 30             for(int j=0; j<M; ++j ){
 31                 for( int k=0; k<M; ++k ){
 32                     C.m[i][j]+=m[i][k]*a.m[k][j];
 33                     C.m[i][j]%=Mod;
 34                 }
 35             }
 36         }
 37         return C;
 38     }
 39     inline Mar operator ^ (int t) const{
 40         Mar B=*this, C;
 41         C.unit( );
 42         while(t){
 43             if(t&1)C=C*B;
 44             B=B*B;
 45             t>>=1;
 46         }
 47         return C;
 48     }
 49 }A;
 50 int a[100005], p[10005],cntp=0;
 51 void getp( )
 52 {
 53     for( int i=2; i<=1e5; ++ i ){
 54         if( !a[i] )p[cntp++]=i;
 55         for( int j=0; j<cntp&&i*p[j]<1e5; ++ j ){
 56             a[i*p[j]]=1;
 57             if( i%p[j]==0 )break;
 58         }
 59     }
 60 }
 61 
 62 int Phi( int x )
 63 {
 64     int res=x;
 65     for( int i=0; i<cntp&&p[i]*p[i]<=x; ++i ){
 66         if( x%p[i]==0 ){
 67             res/=p[i]; res*=(p[i]-1);
 68             while(x%p[i]==0){
 69                 x/=p[i];
 70             }
 71         }
 72     }
 73     if(x>1){
 74         res/=x;res*=(x-1);
 75     }
 76     return res%Mod;
 77 }
 78 int P_M(int a, int b)
 79 {
 80     int res=1;
 81     while(b){
 82         if(b&1)res*=a, res%=Mod;
 83         a*=a, a%=Mod;
 84         b>>=1;
 85     }
 86     return res;
 87 
 88 }
 89 int work( int k )
 90 {
 91     Mar C=A^k;
 92     int res=0;
 93     for( int i=0; i<M; ++i )
 94         res+=C.m[i][i];
 95     return res;
 96 }
 97 int polya(  )
 98 {
 99     int ans=0;
100     for( int i=1; i*i<=N; ++ i ){
101         if( N%i==0 ){
102             if(i*i==N){
103                 ans+=Phi(i)*work(i);
104                 ans%=Mod;
105             }
106             else{
107                 ans+=Phi(i)*work(N/i);
108                 ans%=Mod;
109                 ans+=Phi(N/i)*work(i);
110                 ans%=Mod;
111             }
112         }
113     }
114     return ans;
115 }
116 int main( )
117 {
118     getp();
119     scanf("%d", &T);
120     while (T--){
121         scanf("%d%d%d", &N, &M, &K);
122         A.one();
123         for( int i=0, x, y; i<K; ++i ){
124             scanf("%d%d", &x, &y);
125             A.m[x-1][y-1]=A.m[y-1][x-1]=0;
126         }
127         int ans=polya();
128         int inv=P_M(N%Mod, Mod-2);
129         ans*=inv;ans%=Mod;
130         printf("%d
", ans);
131     }
132     return 0;
133 }

View Code