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后的对角线上的元素和即为所求。

#include <iostream>
#include <cmath>
#include <cstdio>
#include <cstring>
using namespace std;
const int Mod=9973;
int N, M, K, T;
struct Mar
{
    int m[15][15];
    inline void zero( ){
        memset(m, 0, sizeof m);
    }
    inline void one( ){
        zero( );
        for( int i=0; i<M; ++i ){
            for( int j=0; j<M; ++j ){
                m[i][j]=1;
            }
        }
    }
    inline void unit( ){
        zero();
        for( int i=0; i<M; ++ i )
            m[i][i]=1;
    }
    inline Mar operator *(const Mar &a) const {
        Mar C;C.zero();
        for( int i=0; i<M; ++i ){
            for(int j=0; j<M; ++j ){
                for( int k=0; k<M; ++k ){
                    C.m[i][j]+=m[i][k]*a.m[k][j];
                    C.m[i][j]%=Mod;
                }
            }
        }
        return C;
    }
    inline Mar operator ^ (int t) const{
        Mar B=*this, C;
        C.unit( );
        while(t){
            if(t&1)C=C*B;
            B=B*B;
            t>>=1;
        }
        return C;
    }
}A;
int a[100005], p[10005],cntp=0;
void getp( )
{
    for( int i=2; i<=1e5; ++ i ){
        if( !a[i] )p[cntp++]=i;
        for( int j=0; j<cntp&&i*p[j]<1e5; ++ j ){
            a[i*p[j]]=1;
            if( i%p[j]==0 )break;
        }
    }
}

int Phi( int x )
{
    int res=x;
    for( int i=0; i<cntp&&p[i]*p[i]<=x; ++i ){
        if( x%p[i]==0 ){
            res/=p[i]; res*=(p[i]-1);
            while(x%p[i]==0){
                x/=p[i];
            }
        }
    }
    if(x>1){
        res/=x;res*=(x-1);
    }
    return res%Mod;
}
int P_M(int a, int b)
{
    int res=1;
    while(b){
        if(b&1)res*=a, res%=Mod;
        a*=a, a%=Mod;
        b>>=1;
    }
    return res;

}
int work( int k )
{
    Mar C=A^k;
    int res=0;
    for( int i=0; i<M; ++i )
        res+=C.m[i][i];
    return res;
}
int polya(  )
{
    int ans=0;
    for( int i=1; i*i<=N; ++ i ){
        if( N%i==0 ){
            if(i*i==N){
                ans+=Phi(i)*work(i);
                ans%=Mod;
            }
            else{
                ans+=Phi(i)*work(N/i);
                ans%=Mod;
                ans+=Phi(N/i)*work(i);
                ans%=Mod;
            }
        }
    }
    return ans;
}
int main( )
{
    getp();
    scanf("%d", &T);
    while (T--){
        scanf("%d%d%d", &N, &M, &K);
        A.one();
        for( int i=0, x, y; i<K; ++i ){
            scanf("%d%d", &x, &y);
            A.m[x-1][y-1]=A.m[y-1][x-1]=0;
        }
        int ans=polya();
        int inv=P_M(N%Mod, Mod-2);
        ans*=inv;ans%=Mod;
        printf("%d
", ans);
    }
    return 0;
}