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;
}