[HDU2865]Birthday Toy（置换群+dp+矩乘）

题目描述

传送门

题解

首先中间放一种颜色，还剩k-1种
设f(i)表示共i个点染k-1中颜色，相邻的两个（首尾）不能相同的方案数
那么f(i)=f(i−1)∗(k−3)+f(i−2)∗(k−2)
也就是说，在f(i−1)最后插入一个不同于首也不同于尾的颜色，在f(i−2)最后插入一个和首相同以及一个和首不同的颜色，这样就可以递推了
然后发现如果循环节数量为1的话是不行的
然后根据定理和一点数论的知识
ans=kn∑i=1n[gcd(i,n)≠1]f(gcd(i,n))
=kn[∑i=1nf(gcd(i,n))−φ(n)∗(k−1)]
=kn[∑d|nφ(nd)f(d)−φ(n)∗(k−1)]
n√枚举约数，n√求φ，然后用矩阵加速dp求f(i)

想了好久好久啊….想了很多不靠谱的方法，n倍经验了…
卡了一波丧心病狂的常数…

代码

#include<algorithm>#include<iostream>#include<cstring>#include<cstdio>#include<cmath>using namespace std;#define LL long long#define N 1005const LL Mod=1000000007LL;int n;LL k,f1,f2,f3,ans;struct data{LL a[4][4];}unit,st,A;int prime[100005],p[100005],phi[100005];void get(){    phi[1]=1;    for (int i=2;i<=100000;++i)    {        if (!p[i]) {prime[++prime[0]]=i;phi[i]=i-1;}        for (int j=1;j<=prime[0]&&i*prime[j]<=100000;++j)        {            p[i*prime[j]]=1;            if (i%prime[j]==0) {phi[i*prime[j]]=phi[i]*prime[j];break;}            else phi[i*prime[j]]=phi[i]*phi[prime[j]];        }    }}LL fast_pow(LL a,int p){    LL ans=1LL;    for (;p;p>>=1,a=a*a%Mod)        if (p&1)            ans=ans*a%Mod;    return ans;}data cheng(data a,data b){    data ans;memset(ans.a,0,sizeof(ans.a));    LL *w;    for (int i=1;i<=3;++i)        for (int j=1;j<=3;++j)        {            w=&(ans.a[i][j]);            (*w)=((*w)+a.a[i][1]*b.a[1][j])%Mod;            (*w)=((*w)+a.a[i][2]*b.a[2][j])%Mod;            (*w)=((*w)+a.a[i][3]*b.a[3][j])%Mod;            }    return ans;}data mapow(data a,int p){    data ans=unit;    for (;p;p>>=1,a=cheng(a,a))        if (p&1)            ans=cheng(ans,a);    return ans;}LL _phi(int x){    if (x<=100000) return (LL)phi[x];    LL ans=x;    for (int i=1;x>1&&i<=prime[0]&&prime[i]*prime[i]<=x;++i)        if (x%prime[i]==0)        {            ans=ans*(prime[i]-1)/prime[i];            while (x%prime[i]==0) x/=prime[i];        }    if (x>1) ans=ans*(x-1)/x;    return ans%Mod;}LL F(int d){    if (d==1) return f1;    if (d==2) return f2;    if (d==3) return f3;    data now=mapow(A,d-3);    now=cheng(st,now);    return now.a[1][1];}int main(){    get();    unit.a[1][1]=unit.a[2][2]=unit.a[3][3]=1LL;    while (~scanf("%d%I64d",&n,&k))    {        ans=0LL;        A.a[1][1]=k-3LL,A.a[2][1]=k-2LL,A.a[1][2]=1LL,A.a[2][3]=1LL;        f1=k-1LL;f2=(k-1L)*(k-2LL)%Mod;f3=(k-1LL)*(k-2LL)%Mod*(k-3LL)%Mod;        st.a[1][1]=f3,st.a[1][2]=f2,st.a[1][3]=f1;        for (int i=1;i*i<=n;++i)            if (n%i==0)            {                ans+=_phi(n/i)*F(i)%Mod;                if (i*i!=n) ans+=_phi(i)*F(n/i)%Mod;                if (ans>=Mod) ans-=Mod;            }        ans=ans-_phi(n)*(k-1LL);ans=(ans%Mod+Mod)%Mod;        ans=ans*k%Mod*fast_pow(n,Mod-2LL)%Mod;        printf("%I64d\n",ans);    }}