HDU 2865 Birthday Toy(Polya+矩阵乘法+dp)

题意：要构成一个长度为n的项链，中间有个大珠子，有k种颜色，问不等价的涂色方案，并且相邻的两个珠子不能同色，通过旋转得到的被认为是等价的。

思路：这题主要的问题就是如何处理相邻珠子不能同色的问题，中间的大珠子好解决，一开始给它分配个颜色就好了，然后用k-1种颜色涂色。思考旋转后置换的效果，为了要求等价，同一循环节颜色要相同，另外还要满足相邻的颜色不能相同，其实这等同于对和循环节数量相等的珠子按这样的要求涂颜色（不明白可以画一画）。那么不妨枚举一下第一个珠子的颜色，dp[i][j]表示第i个珠子是（否）是给定的颜色的涂色方案数，0表示不同色，1表示同色，那么就可以得到dp[i][0]=dp[i-1][0]*(k-3)+dp[i-1][1]*(k-2),dp[i][1]=dp[i-1][0],这样就得到了递推式，用矩阵乘法可以快速求出，最终合法方案就是最后一个珠子的颜色等于给定颜色的方案数，即dp[n][1]，由于第一个珠子颜色用哪个结果都相同，则结果乘k-1就行。

代码：

#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=1000000007;const int maxn=40000+10;//矩阵乘法相关struct Matrix{    ll mat[2][2];    void Init()    {        mat[0][0]=mat[1][1]=1;        mat[1][0]=mat[0][1]=0;    }    void clear() {memset(mat,0,sizeof(mat));}}mx;Matrix operator *(const Matrix &a,const Matrix &b){    Matrix c;c.clear();    for(int k=0;k<2;++k)        for(int i=0;i<2;++i)            for(int j=0;j<2;++j)                c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;    return c;}Matrix pow_mat(Matrix x,ll n){    Matrix res;res.Init();    while(n)    {        if(n&1) res=res*x;        x=x*x;        n>>=1;    }    return res;}//数论相关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;        }    }}ll euler_phi(ll n){    ll 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;}ll pow_mod(ll x,ll n){    ll res=1;    while(n)    {        if(n&1) res=res*x%mod;        x=x*x%mod;        n>>=1;    }    return res;}ll inv(ll a){    return pow_mod(a,mod-2);}int n,k;ll ans;int factor[110],num[110],tot;void dfs(int x,int now){    if(x==tot)    {        if(now==1) return ;        Matrix y=pow_mat(mx,now-1);        ans+=y.mat[1][0]*euler_phi(n/now);        ans%=mod;        return ;    }    int tmp=1;    for(int i=0;i<=num[x];++i)    {        dfs(x+1,now*tmp);        tmp*=factor[x];    }}void solve(){    int N=n;    tot=0;    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;}    mx.clear();    mx.mat[0][0]=k-3;mx.mat[1][0]=k-2;    mx.mat[0][1]=1;mx.mat[1][1]=0;    dfs(0,1);    ans=ans*(k-1)%mod;    ans=ans*k%mod;    ans=ans*inv(n)%mod;}int main(){    //freopen("in.txt","r",stdin);    //freopen("out.txt","w",stdout);    getprimes();    while(~scanf("%d%d",&n,&k))    {        ans=0;        solve();        printf("%I64d\n",ans);    }    return 0;}