首页 > 代码库 > HDU 2865 Birthday Toy [Polya 矩阵乘法]

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

 

HDU 2865 Birthday Toy [Polya 矩阵乘法]