[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={k\over n}\sum\limits_{i=1}^n[gcd(i,n)\neq 1]f(gcd(i,n))$
$={k\over n}[\sum\limits_{i=1}^nf(gcd(i,n))-\varphi(n)*(k-1)]$
$={k\over n}{[\sum\limits_{d|n}\varphi({n\over d})f(d)-\varphi(n)*(k-1)]}$
$\sqrt n$枚举约数，$\sqrt n$求$\varphi$，然后用矩阵加速dp求$f(i)$

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

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define LL long long
#define N 1005

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