POJ 2888 Magic Bracelet (Polya+欧拉函数+矩阵乘法)

题意：用m种颜色的珠子串成长度为n的项链，问有多少种不等价的方案，给出的珠子中，有一些珠子是不能同时出现在项链中的。

思路：刚开始看的时候完全没有头绪，主要是给出的限制条件不好处理，后来看了下题解才知道可以用矩阵来搞，把串成的链看成一个路径的话，那么就是从某个颜色的点出发，走n歩然后再回到n点的路径数，这个正好可以用矩阵乘法解决。。。然后对于每个循环节，由于需要涂同样的颜色，可以看成一个珠子，然后求相应长度的方案数就行了。对于求循环节的方法，要枚举因子然后用欧拉函数进行优化。

代码:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<map>
#include<queue>
#include<stack>
#include<cmath>
#include<vector>
#define inf 0x3f3f3f3f
#define Inf 0x3FFFFFFFFFFFFFFFLL
#define eps 1e-9
#define pi acos(-1.0)
using namespace std;
typedef long long ll;
const int mod=9973;
const int maxn=40000+10;
struct Matrix
{
    int mat[11][11],size;
    void Init(int size)
    {
        this->size=size;
        for(int i=0;i<size;++i)
            for(int j=0;j<size;++j)
                mat[i][j]=(i==j);
    }
    void clear() {memset(mat,0,sizeof(mat));}
}mx;
Matrix operator *(const Matrix &a,const Matrix &b)
{
    Matrix c;c.clear();
    c.size=a.size;
    for(int k=0;k<c.size;++k)
        for(int i=0;i<c.size;++i)
            for(int j=0;j<c.size;++j)
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
    return c;
}
Matrix pow_mat(Matrix a,int n)
{
    Matrix res;
    res.Init(a.size);
    while(n)
    {
        if(n&1) res=res*a;
        a=a*a;
        n>>=1;
    }
    return res;
}
int pow_mod(int x,int n)
{
    int res=1;
    x%=mod;
    while(n)
    {
        if(n&1) res=res*x%mod;
        x=x*x%mod;
        n>>=1;
    }
    return res;
}
int inv(int a)
{
    return pow_mod(a,mod-2)%mod;
}
int primes[maxn],pcnt;
bool flag[maxn];
void getprimes()
{
    memset(flag,0,sizeof(flag));
    pcnt=0;
    for(int i=2;i<maxn;++i)
    {
        if(!flag[i])
        {
            primes[pcnt++]=i;
            for(int j=i*i;j<maxn;j+=i)
                flag[j]=true;
        }
    }
}
int euler_phi(int n)
{
    int ans=n;
    for(int i=0;primes[i]*primes[i]<=n;++i)
    {
        if(n%primes[i]==0)
        {
            ans=ans-ans/primes[i];
            while(n%primes[i]==0) n/=primes[i];
        }
    }
    if(n>1) ans=ans-ans/n;
    return ans%mod;
}
int factor[110],num[110],tot,ans;
void dfs(int x,int now,int n)
{
    if(x==tot)
    {
        int val=0;
        Matrix y=pow_mat(mx,now);
        for(int i=0;i<mx.size;++i)
            val=(val+y.mat[i][i])%mod;
        ans=(ans+euler_phi(n/now)*val)%mod;
        return ;
    }
    int tmp=1;
    for(int i=0;i<=num[x];++i)
    {
        dfs(x+1,now*tmp,n);
        tmp*=factor[x];
    }
}
void solve(int n)
{
    ans=0;tot=0;
    int N=n;
    for(int i=0;primes[i]*primes[i]<=N;++i)
    {
        if(N%primes[i]==0)
        {
            factor[tot]=primes[i];
            num[tot]=0;
            while(N%primes[i]==0) {num[tot]++;N/=primes[i];}
            tot++;
        }
    }
    if(N>1) {factor[tot]=N;num[tot++]=1;}
    dfs(0,1,n);
    ans=ans*inv(n)%mod;
}
int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    getprimes();
    int t,n,m,k;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d%d",&n,&m,&k);
        mx.size=m;
        mx.clear();
        for(int i=0;i<mx.size;++i)
            for(int j=0;j<mx.size;++j)
                mx.mat[i][j]=1;
        int u,v;
        while(k--)
        {
            scanf("%d%d",&u,&v);
            mx.mat[u-1][v-1]=0;
            mx.mat[v-1][u-1]=0;
        }
        solve(n);
        printf("%d\n",ans);
    }
    return 0;
}