POJ 2888 Magic Bracelet （经典Polya+矩阵）

转载请注明出处，谢谢http://blog.csdn.net/ACM_cxlove?viewmode=contents by---cxlove

题目：同样是S个物品的环，C种颜色，旋转视为相同。但是有一些颜色限制，规定某种颜色不能相邻。

http://poj.org/problem?id=2888

难点在于颜色的限制。 巧妙的运用矩阵，a->b->c->a，表示一个长度为3的一种方案，而且最终是回到起点的。那么便是从一个点出发，经过L步，返回起始点，经典的矩阵问题，通过颜色限制求出矩阵，由矩阵快速幂乘，枚举起点颜色，就能得到一定长度的所有方案数，也就是矩阵的对角线和。

最后除以置换群数，记得取逆元。HDU 2865有着类似的应用，不过范围更大，不能直接用矩阵，必须推出公式

#include<iostream>
#include<cstring>
#include<queue>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define N 1000000000
#define inf 1<<29
#define MOD 9973
#define LL long long
using namespace std;
struct Matrix{
    int m[15][15];
}init;
int s,c,k;
bool flag[40000]={0};
int prime[40000],cnt=0;
Matrix operator*(Matrix m1,Matrix m2){
    Matrix ans;
    for(int i=0;i<c;i++)
        for(int j=0;j<c;j++){
            ans.m[i][j]=0;
            for(int k=0;k<c;k++)
                ans.m[i][j]=(ans.m[i][j]+m1.m[i][k]*m2.m[k][j])%MOD;
        }
    return ans;
}
Matrix operator^(Matrix m1,int b){
    Matrix ans;
    for(int i=0;i<c;i++)
        for(int j=0;j<c;j++)
            ans.m[i][j]=(i==j);
    while(b){
        if(b&1)
            ans=ans*m1;
        m1=m1*m1;
        b>>=1;
    }
    return ans;
}
void Prime(){
    for(int i=2;i<=sqrt(N+1.0);i++){
        if(flag[i]) continue;
        prime[cnt++]=i;
        for(int j=2;j*i<=sqrt(N+1.0);j++)
            flag[i*j]=true;
    }
}
int PowMod(int a,int b){
    a%=MOD;
    int ret=1;
    while(b){
        if(b&1)  ret=(ret*a)%MOD;
        a=(a*a)%MOD;
        b>>=1;
    }
    return ret;
}
int Eular(int n){
    int ret=1;
    for(int i=0;i<cnt&&prime[i]*prime[i]<=n;i++){
        if(n%prime[i]==0){
            n/=prime[i];ret*=prime[i]-1;
            while(n%prime[i]==0){n/=prime[i];ret=(ret*prime[i])%MOD;}
        }
    }
    if(n>1) ret*=n-1;
    return ret%MOD;
}
void debug(Matrix t){
    for(int i=0;i<c;i++){
        for(int j=0;j<c-1;j++)
            printf("%d ",t.m[i][j]);
        printf("%d\n",t.m[i][c-1]);
    }
}
int slove(int n){
    Matrix temp=init^n;
    int ans=0;
    for(int i=0;i<c;i++)
        ans=(ans+temp.m[i][i])%MOD;
    return ans;
}
int Polya(){
    int i,ans=0;
    for(i=1;i*i<s;i++)
        if(s%i==0)
           ans=(ans+Eular(i)*slove(s/i)+Eular(s/i)*slove(i))%MOD;
    if(i*i==s) ans=(ans+Eular(i)*slove(i))%MOD;
    return (ans*PowMod(s%MOD,MOD-2))%MOD;
}
int main(){
    int t;
    Prime();
    scanf("%d",&t);
    while(t--){
        scanf("%d%d%d",&s,&c,&k);
        for(int i=0;i<c;i++)
            for(int j=0;j<c;j++)
                init.m[i][j]=1;
        while(k--){
            int a,b;
            scanf("%d%d",&a,&b);
            a--;b--;
            init.m[a][b]=init.m[b][a]=0;
        }
        printf("%d\n",Polya());
    }
    return 0;
}