POJ 2888 Magic Bracelet（polya+矩阵快速幂）

Description
用m种颜色给一串长度为n的项链染色，其中有k对颜色不能相邻，旋转相同看作同一种方案，问方案数，结果模9973
Input
第一行为一整数T表示用例组数，每组用例第一行为三个整数n,m,k分别表示颜色数，珠子数以及限制数，之后k行每行两个整数a和b表示颜色a和颜色b不能相邻
Output
对于每组用例，输出染色方案数，结果模9973
Sample Input
4
3 2 0
3 2 1
1 2
3 2 2
1 1
1 2
3 2 3
1 1
1 2
2 2
Sample Output
4
2
1
0
Solution
polya，置换数和每种置换对应的轮换数和没有限制时相同，问题在于对于每一种置换有多少种染色方案，设轮换数为t，m种颜色的关系矩阵为M，令A=M^t，那么第i种颜色经过A[i][i]种可行路径就会变成自己，所以总方案数为tr(A)
Code

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define maxm 55555
#define maxn 11
#define mod 9973
int prime[maxm],is_prime[maxm],res;
void get_prime(int n)
{
    res=0;
    memset(is_prime,0,sizeof(is_prime));
    for(int i=2;i<n;i++)
        if(!is_prime[i])
        {
            prime[res++]=i;
            for(int j=i;j<n;j+=i)
                is_prime[j]=1;
        }
}
int mod_pow(int a,int b,int p)
{
    int ans=1;
    a%=p;
    while(b)
    {
        if(b&1) ans=(ans*a)%p;
        a=(a*a)%p;
        b>>=1;
    } 
    return ans;
}
int get_euler(int n)
{
    int ans=n;
    for(int i=0;i<res&&prime[i]*prime[i]<=n;i++)
        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;
} 
struct Mat
{
    int mat[maxn][maxn];
    int row,col;
};
Mat mod_mul(Mat a,Mat b,int p)
{
    Mat ans;
    ans.row=a.row;
    ans.col=b.col;
    memset(ans.mat,0,sizeof(ans.mat));
    for(int i=0;i<ans.row;i++)      
        for(int k=0;k<a.col;k++)
            if(a.mat[i][k])
                for(int j=0;j<ans.col;j++)
                {
                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                    ans.mat[i][j]%=p;
                }
    return ans;
}
Mat mod_pow(Mat a,int k,int p) 
{
    Mat ans;
    ans.row=a.row;
    ans.col=a.col;
    for(int i=0;i<a.row;i++)
        for(int j=0;j<a.col;j++)
            ans.mat[i][j]=(i==j);
    while(k)
    {
        if(k&1)ans=mod_mul(ans,a,p);
        a=mod_mul(a,a,p);
        k>>=1;
    }
    return ans;
}
int get_track(Mat M,int n,int p)
{
    M=mod_pow(M,n,p);
    int ans=0;
    for(int i=0;i<M.row;i++)
        ans=(ans+M.mat[i][i])%p;
    return ans;
}
int main()  
{  
    get_prime(maxm);
    int T,n,m,k;   
    scanf("%d",&T);  
    while(T--)  
    {  
        scanf("%d%d%d",&n,&m,&k);
        Mat M;
        M.row=M.col=m;
        for(int i=0;i<m;i++)
            for(int j=0;j<m;j++)
                M.mat[i][j]=1;
        while(k--)
        {
            int a,b;
            scanf("%d%d",&a,&b);
            M.mat[a-1][b-1]=M.mat[b-1][a-1]=0;
        }
        int ans=0;
        for(int i=1;i*i<=n;i++)
            if(n%i==0)
            {
                if(i*i!=n)
                    ans=(ans+get_euler(i)%mod*get_track(M,n/i,mod)%mod+get_euler(n/i)%mod*get_track(M,i,mod)%mod)%mod;
                else 
                    ans=(ans+get_euler(i)%mod*get_track(M,n/i,mod)%mod)%mod;
            }
        ans=ans*mod_pow(n,mod-2,mod)%mod;
        printf("%d\n",ans);
    }  
    return 0;  
}