POJ2888 Magic Bracelet【Burnside引理 + DP(矩阵加速)求不动点】

题目描述：

大小为n，可用颜色数为m的一串珠子，其中有k对限制某两种颜色不能相邻。问在忽略中心旋转产生的重复下的不同的染色方案数。
$n\le 10^9,m\le 10,k\le \frac{m(m-1)}{2}$

题目分析：

如果没有限制就是裸的Polya，直接phi乘上颜色种数的几次幂求和。

但是有限制就会导致在循环$i$步的置换下形成的长度为$gcd(i,n)$的循环节它的不动点数不再是$m^{gcd(i,n)}$。

考虑DP求$f[i][j]$表示考虑了前$i$个点，第$i$个点颜色为$j$的方案数，容易写出转移并用矩阵加速。但是注意第1个点和第$gcd$个点也是相邻的，我们可以先钦定第0个位置为颜色$c$，最后答案就加上$f[n][c]$。
这样看起来需要做$m$次矩阵快速幂，但是由于矩阵乘法的左矩阵的一行不会对另一行做贡献，所以我们可以直接令初始矩阵的$a[i][i]=1$，最后加上每一行的$a[i][i]$即可，相当于把多次行矩阵的加速合在一起了。

Code：

#include<cstdio>
#include<cstring>
const int mod = 9973;
int n,m,k,ans,d[3005],phi[3005],cnt;
struct Mat{
    int s[10][10];
    Mat(){memset(s,0,sizeof s);}
    void unit(){memset(s,0,sizeof s);for(int i=0;i<m;i++) s[i][i]=1;}
    Mat operator * (const Mat &B)const{
        Mat ret;
        for(int k=0;k<m;k++)
            for(int i=0;i<m;i++) if(s[i][k])
                for(int j=0;j<m;j++) if(B.s[k][j])
                    ret.s[i][j]+=s[i][k]*B.s[k][j];
        for(int i=0;i<m;i++) for(int j=0;j<m;j++) ret.s[i][j]%=mod;
        return ret;
    }
}f,g,G;
inline void Add(int x,int t){
    int now=cnt;
    for(int i=1;i<=now;i++){
        d[++cnt]=d[i]*x,phi[cnt]=phi[i]*(x-1);
        for(int k=2;k<=t;k++) d[cnt+1]=d[cnt]*x,phi[cnt+1]=phi[cnt]*x,cnt++;
    }
}
inline int MatPow(int N){
    f.unit(),g=G;
    for(;N;N>>=1,g=g*g) if(N&1) f=f*g;
    int ret=0;
    for(int i=0;i<m;i++) ret+=f.s[i][i];
    return ret%mod;
}
inline int Pow(int a,int b){int s=1;for(;b;b>>=1,a=a*a%mod) if(b&1) s=s*a%mod;return s;}
int main()
{
    scanf("%*d");
    while(~scanf("%d%d%d",&n,&m,&k)){
        ans=0,d[cnt=1]=1,phi[1]=1;
        int x=n,y,t;
        for(int i=2;i*i<=x;i++) if(x%i==0){t=0;while(x%i==0) x/=i,t++;Add(i,t);}
        if(x>1) Add(x,1);
        for(int i=0;i<m;i++) for(int j=0;j<m;j++) g.s[i][j]=1;
        for(int i=1;i<=k;i++) scanf("%d%d",&x,&y),x--,y--,g.s[x][y]=g.s[y][x]=0;
        G=g;
        for(int i=1;i<=cnt;i++) ans=(ans+1ll*phi[i]*MatPow(n/d[i]))%mod;
        printf("%d\n",ans*Pow(n%mod,mod-2)%mod);
    }
}