HDU 2865 Birthday Toy [Polya 矩阵乘法]

传送门

题意：

相邻珠子不能相同，旋转等价。$n$个珠子$k$中颜色，求方案数

---

 

首先中间珠子$k$种选择，$k--$
如果没有相邻不同的限制，就和$POJ\ 2154$一样了
$|C(f)|=k^{\#(f)}$
但是有了相邻不同的限制，每种循环的颜色就不能任意选择了
旋转等价循环个数是$gcd(n,i)$，同一个循环的元素相差$i$步
容易得到只要满足长度$gcd(n,i)$的一段相邻颜色不同整个环就不同了，因为这样的一段正好每个循环有一个元素
考虑$DP$，$f[i]$表示$i$个元素组成的环染色方案数
$f[i]=(k-2)*f[i-1]+(k-1)*f[i-2]$
因为这时候$i-1$是可以和$1$相同的，那样$i$有$k-1$种选择，所以加上后面的一块
显然可以用矩阵快速幂
计算的时候使用和和$POJ\ 2154$同样的技巧
最后的式子为：
$\frac{k}{n}\sum\limits_{d \mid n}{f(d)*\phi(\frac{n}{d})},\ d \neq 1$

 

注意：$Candy?$把矩阵的构造函数里加上了每个矩阵都初始化为单位矩阵，认为这样就不用在做矩阵快速幂前初始化了。

然后就被坑惨了......矩阵乘法里还需要零矩阵啊啊啊啊啊啊啊 

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N=1e5+5,P=1e9+7;
typedef long long ll;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}
    return x*f;
}
int n;
ll k;
int p[N];
bool notp[N];
void sieve(int n){
    for(int i=2;i<=n;i++){
        if(!notp[i]) p[++p[0]]=i;
        for(int j=1;j<=p[0]&&i*p[j]<=n;j++){
            notp[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
}
inline int phi(int n){
    int re=n,m=sqrt(n);
    for(int i=1;i<=p[0]&&p[i]<=m&&p[i]<=n;i++) if(n%p[i]==0){
        re=re/p[i]*(p[i]-1);
        while(n%p[i]==0) n/=p[i];
    }
    if(n>1) re=re/n*(n-1);
    return re%P;
}
struct Matrix{
    ll a[2][2];
    ll* operator [](int x){return a[x];}
    Matrix(){a[0][0]=a[1][1]=a[0][1]=a[1][0]=0;}
    void ini(){a[0][0]=a[1][1]=1;}
}a,b;
Matrix operator *(Matrix a,Matrix b){
    Matrix c;
    for(int k=0;k<2;k++)
        for(int i=0;i<2;i++) if(a[i][k])
            for(int j=0;j<2;j++) if(b[k][j])
                (c[i][j]+=a[i][k]*b[k][j])%=P;
    return c;
}
Matrix operator ^(Matrix a,int b){
    Matrix re;re.ini();
    for(;b;b>>=1,a=a*a)
        if(b&1) re=re*a;
    return re;
}
ll F[5];
ll f(int x){
    if(x<=3) return F[x];
    Matrix re=a^(x-3);
    re=re*b;
    return re[0][0];
}
inline void mod(int &x){if(x>=P) x-=P;}
inline ll Pow(ll a,int b){
    ll re=1;
    for(;b;b>>=1,a=a*a%P)
        if(b&1) re=re*a%P;
    return re;
}
inline ll Inv(ll a){return Pow(a,P-2);}
void solve(){
    int m=sqrt(n),ans=0;
    for(int i=1;i<=m;i++) if(n%i==0){
        if(i!=1) mod(ans+= f(i)*phi(n/i)%P);
        if(i*i!=n) mod(ans+= f(n/i)*phi(i)%P);
    }
    printf("%lld\n",ans*Inv(n)%P*(k+1)%P);
}
int main(){
    freopen("in","r",stdin);
    sieve(32000);
    while(scanf("%d%lld",&n,&k)!=EOF){
        k--;
        F[1]=k;F[2]=k*(k-1)%P;F[3]=k*(k-1)%P*(k-2)%P;
        a[0][0]=k-2; a[0][1]=k-1;
        a[1][0]=1;     a[1][1]=0;
        b[0][0]=F[3];b[0][1]=0;
        b[1][0]=F[2];b[1][1]=0;
        solve();
    }
}