[HDU2481]Toy（置换+矩阵乘法）

=== ===

这里放传送门

=== ===

题解

一眼看到这个题的时候想到了BZOJ1002轮状病毒。。。实际上做这个题的时候也用到了那个题的结论。

结论：设f(n)为不同的（不考虑旋转同构）的n轮状病毒的个数，则f(n)=3∗f(n−1)−f(n−2)+2

这玩意儿当时是用基尔霍夫矩阵推出来的。。
但是这个题的话要考虑旋转同构诶，如何和那个题的结论联系起来呢？

Burnside引理：设有置换群G={a1,a2...ag}，c(ai)为置换ai作用下不动点的个数，那么这个置换群作用下等价类的个数为1G∑i=1gc(ai)

那么关键就是如何求出在每种置换作用下不动点的个数，也就是经过这个置换作用以后不变的方案数目。设旋转i次的置换为ai，那么显然置换总数为n，ai的循环节个数为gcd(i,n)。
可以发现如果这个图被分成了i个循环节，那么属于第i个循环节的点一定和属于第i%n+1和第i%n−1个循环节（啊这里什么减出负数来的细节就不管了反正就是要表达一个循环中的前一个点和后一个点的意思啦）的点是挨着的。

那么如果要构造在置换ai作用下的同构方案的话，显然如果某个属于循环节x的点和某个属于循环节y之间的点有边相连，那么所有属于循环节x的点都必须和对应位置循环节y的点相连，不然旋转一下的话看起来就不一样了。
这说明这个连边方案只和循环节的个数有关，属于同一个循环节的点可以看成一个整体，而如果把循环节看成整体的话就相当于可以随便连边不用考虑什么同构限制了。那这就跟轮状病毒那个题是一样了。
设n轮状病毒的方案数为f(n)，那么我们要求的结果就是∑i=1nf(gcd(i,n))。
用数学方法化一下式子就可以化成∑d|nφ(d)f(nd)，用O(n√)的时间就可以求解。

注意的事情是这里的模数不一定和n互质，在做除法的时候没法求逆元，那只能先把m乘上n，最后除以n再模m，而这就造成了模数非常大，直接乘可能炸long long，所以要用快速乘。

代码

#include<cstdio>#include<cstring>#include<algorithm>using namespace std;int n,prm[1000010];long long ans,m,f[5],Mod,phi[1000010];bool ext[1000010];inline long long mul(long long a,long long t){    long long ans=0;    if (a<0) a+=Mod;    if (t<0) t+=Mod;    while (t!=0){        if (t&1) ans=ans+a;        if (ans>=Mod) ans-=Mod;        a=a+a;t>>=1;        if (a>=Mod) a-=Mod;    }    return ans;}struct hxy{    long long s[5][5];    hxy(){memset(s,0,sizeof(s));}    void clear(){        memset(s,0,sizeof(s));        for (int i=1;i<=3;i++) s[i][i]=1;    }    void get(){        memset(s,0,sizeof(s));        s[1][1]=1;s[3][2]=1;s[1][3]=2;        s[2][3]=-1;s[3][3]=3;    }    hxy operator * (const hxy &a){        hxy c;        long long *w;        for (int i=1;i<=3;i++)          for (int j=1;j<=3;j++){              w=&(c.s[i][j]);              (*w)=((*w)+mul(s[i][1],a.s[1][j]))%Mod;              (*w)=((*w)+mul(s[i][2],a.s[2][j]))%Mod;              (*w)=((*w)+mul(s[i][3],a.s[3][j]))%Mod;          }        return c;    }}w;hxy powww(hxy a,int t){    hxy ans;    ans.clear();    while (t!=0){        if (t&1) ans=ans*a;        a=a*a;t>>=1;    }    return ans;}void get_prime(int N){    phi[1]=1;    for (int i=2;i<=N;i++){        if (ext[i]==false){            prm[++prm[0]]=i;            phi[i]=i-1;        }        for (int j=1;j<=prm[0];j++){            if ((long long)i*prm[j]>N) break;            ext[i*prm[j]]=true;            if (i%prm[j]==0){                phi[i*prm[j]]=phi[i]*prm[j];break;            }else phi[i*prm[j]]=phi[i]*phi[prm[j]];        }    }}long long PHI(int N){    if (N<=1000000) return phi[N]%Mod;    long long ans=N;    for (int i=1;i<=prm[0];i++){        int v=prm[i];        if ((long long)v*v>N) break;        if (N%v==0) ans=ans-ans/v;        while (N%v==0) N/=v;        if (N==1) break;    }    if (N!=1) ans=ans-ans/N;    return ans%Mod;}long long F(int N){    long long h[5];    if (N<=2) return f[N+1];    memset(h,0,sizeof(h));    w.get();    w=powww(w,N-2);    for (int i=1;i<=3;i++)      for (int j=1;j<=3;j++)        h[i]=(h[i]+mul(f[j],w.s[j][i]))%Mod;    return h[3];}int main(){    get_prime(1000000);    while (scanf("%d%d",&n,&m)!=EOF){        ans=0;Mod=(long long)m*n;        f[1]=1;f[2]=1;f[3]=5;        for (int i=1;i*i<=n;i++)          if (n%i==0)            if (i*i==n) ans=(ans+mul(PHI(i),F(i)))%Mod;            else{                ans=(ans+mul(PHI(i),F(n/i)))%Mod;                ans=(ans+mul(PHI(n/i),F(i)))%Mod;            }        ans=(ans+Mod)%Mod;        ans=(ans/n)%m;        printf("%I64d\n",ans);    }    return 0;}