poj 2888 Magic Bracelet

一道很不错的题目，这里加入连接限制同时还考察优化，优化方法同上。连接限制如何处理？注意到项链个数很少，因此可以建图，然后分别求出每种 颜色连接n个珠子后回到自身的方案数，累加即可，这里可以用矩阵快速幂求解。

 

#include<iostream>
 #include<cstdio>
 #include<cstdlib>
 #include<algorithm>
 #include<cmath>
 #include<queue>
 #include<set>
 #include<map>
 #include<cstring>
 #include<vector>
 #include<string>
 #define LL long long
 #define Mod 9973
 using namespace std;
 int prime[4500],cnt;
 class Matrix
 {
 public:
       int matrix[15][15];
       void Zero(  ){ memset( matrix ,0 , sizeof( matrix ) ); }
       int n;
       Matrix operator *( Matrix &b );    
 };
 void Prime(  )
 {
     bool hash[20024] = {0};
     for( int i = 3 ; i <= 200 ; i += 2 ){
             if( !hash[i>>1] ){
                 int x = i << 1;
                 for( int j = i*i; j <= 40000; j+=x )
                      hash[j>>1] = true;
                 }
        }    
     cnt = 0;
     prime[cnt++] = 2;
     for( int i = 1; i <= 20000; i ++ ){
         if( !hash[i] ) prime[cnt++] = i <<1|1;
         }
 }
 Matrix Matrix::operator *( Matrix &b )
 {
       Matrix t;
       t.Zero();
       for( int i = 0 ;i < n; i ++ ){
            for( int k = 0 ; k < n ; k ++ ){
                 if( matrix[i][k] ){
                     for( int j = 0 ; j < n ; j ++ ){
                          if( b.matrix[k][j] )
                              t.matrix[i][j] = ( t.matrix[i][j] + matrix[i][k]*b.matrix[k][j] )%Mod;
                     }    
                 }        
            }    
       }    
       t.n = n;    
       return t;
 }
 int Eular( int n ){
     int ans = n;
     for( int i = 0 ; i < cnt ; i ++ ){
         if( n < prime[i] ) break;
         if( n % prime[i] == 0 ){
             ans = ans/prime[i]*(prime[i]-1);
             while( n % prime[i] == 0 ) {
                 n /= prime[i];
                 }
             }
         }
         if( n != 1 ) ans = ans/n*(n-1);
         return ans%Mod;
 }
 int Pow( int c, int k ){
     int ans = 1;
     c %= Mod;
     while( k ){
         if( k&1 ) ans = (ans*c)%Mod;
         c = (c*c)%Mod ;
         k >>= 1;
         }
     return ans;
     }
 int Solve( Matrix A ,int n  ){
     Matrix t;
     t.Zero();
     for( int i = 0; i < A.n; i ++ ) t.matrix[i][i] = 1;
     t.n = A.n;
     while( n ){
         if( n & 1 ) t =t*A;
         A = A*A;
         n >>= 1;
         }
     int sum =0;
     for( int i = 0 ; i < A.n ; i ++ )
          sum = (sum +t.matrix[i][i])%Mod;
     return sum;
     }
 int Polya( Matrix A,int n ,int k ){
     int t = (int)sqrt( (double)n ),sum  = 0;
     for( int i = 1; i <= t; i ++  ){
         if( n % i == 0 ){
     //        printf( "%d %d %d %d %d\n",sum,Eular(i),Eular(n/i ),Solve( A,i ),Solve( A , n /i ) );
             if( i*i == n )
                sum = ( sum + Eular( i )*Solve(A, i ) )%Mod;
              else
                  sum = ( sum + Eular( i )*Solve( A,n/i )+ Eular( n/i )*Solve(A , i ) )%Mod;
             }
         }
     sum = (sum * Pow( n , Mod - 2 ))%Mod;
     return sum;
     }
 int main(  )
 {
     int n , m , T,k,a,b;
     Prime();
     int p = 
     scanf( "%d",&T );
     while( T-- ){
           scanf( "%d %d %d",&n,&m,&k );
           Matrix A;A.Zero();
           A.n = m;
           for( int i = 0 ; i < m ; i ++ ){
             for( int j = 0 ; j < m ; j ++ ){
                 A.matrix[i][j] = 1;
                 }
             }
         for( int i = 0 ; i < k ;i ++ ){
              scanf( "%d %d",&a,&b );
              A.matrix[a-1][b-1] = 0;
              A.matrix[b-1][a-1] = 0;
         }
         printf( "%d\n",Polya( A,n ,k ) );
     }  
     //system( "pause" );
     return 0;
 }