【hdu2481】Toy，burnside引理+矩阵乘法

传送门
思路：
快把我做哭了TAT
从昨天上午开始想，搞了一下午有一个点没有想明白
TA爷看过题后想了5min貌似就爆正解了TAT
下面我就来讲一讲~
一开始先想没有置换情况下的方案数
手玩无果后打了表……
然后就是
$F(1)=1,F(2)=5,F(n)=3×F(n-1)-F(n-2)+2(n≥3)$
mrazer：这不就是“轮状病毒”吗
咦好像真的是啊……
当时那道题我好像也是打表搞得啊……
说说比较好理解的一种证明
设$F(n)$是带环的方案数，即n号点可以连到1号点
$f(n)$是不带环的方案数，即n号点不能连到1号点
$f(n)= \begin{cases} 1&(n=0,1)\\ 2×f(n-1)+\sum^n_{i=2}f(n-i)&(n≥2)\\ \end{cases}$
$F(n)=\sum^n_{i=1}i^2×f(n-i)$
先来解释一下$F(n)$
考虑把图形外面的n条边想象成分开的一段段
枚举第一个点的段的大小i(1~n)，因为这个段是独立的，想在整体中就必须向中心点连边，所以再枚举这个段向中心点连边的点，数量也是i
剩下的话就是递归问题了，不过要求剩下的段不能成环，就是$f(n-i)$了
$f(n)$的话，考虑前n-1个点已经排好了，加入第n个点使其仍然是树，有三种选择，把n直接连在n-1所在的段上，把n直接与中心点相连，把n连在n-1所在的段上并与中心点相连（也就意味着n-1所在的段形态确定，除n以外不与中心点相连），三种选择的方案分别是$f(n-1)$,$f(n-1)$,$\sum^n_{i=2}f(n-i)$，最后一种相当于枚举n-1所在段的大小i(1-n)，方案数就取决于剩下的n-1-i个点了
但我们发现上面的式子求解在$n≤10^9$下不适用，所以要化简
把$f(n)\ \ (n≥3)$化简一下就成了这个样子

$$f(n)=f(n-1)+\sum^{n-1}_{i=2}f(i)\ \ \ \ \ \ ①$$
$$f(n-1)=f(n-2)+\sum^{n-2}_{i=2}f(i)\ \ \ \ \ ②$$
①-②,移项得
$f(n)=3f(n-1)-f(n-2) \ \ (n≥3)$
（以下过程我写的有些鬼畜， 推荐由结论推出条件）
$$∵\sum^n_{i=n-2}i^2×f(n-i)-3(n-1)^2-2(n-2)^2-2=0$$
$$∴F(n)=\sum^n_{i=1}i^2×f(n-i)-\sum^n_{i=n-2}i^2×f(n-i)+3(n-1)^2+2(n-2)^2+2$$
$$=\sum^{n-3}_{i=1}i^2×f(n-i)+3(n-1)^2+2(n-2)^2+2$$
$∵n-i≥3\\∴f(n-i)=3(n-i-1)-f(n-i-2)$
$$=\sum^{n-3}_{i=1}i^2×[3f(n-i-1)-f(n-i-2)]+3(n-1)^2×f(0)+(n-2)^2×[3f(2)-f(1)]+2$$
$$=\sum^{n-2}_{i=1}i^2×[3f(n-i-1)-f(n-i-2)]+3(n-1)^2×f(0)+2$$
$$=3\sum^{n-1}_{i=1}i^2×f(n-1-i)-\sum^{n-2}_{i=1}i^ 2×f(n-2-i)+2$$
$$=3F(n-1)-F(n-2)+2$$
终于搞出了上面的结论了……
对于这种线性递推式就可以矩乘了……
初始矩阵: $\begin{bmatrix} F[2]&F[1]&2\\ 0&0&0\\ 0&0&0\\ \end{bmatrix}$
转移矩阵: $\begin{bmatrix} 3&1&0\\ -1&0&0\\ 1&0&1\\ \end{bmatrix}$
也有是用行列式、基尔霍夫矩阵直接秒的……蒟蒻不太会，准备找个时间去学习一下O__O
然后就是套burnside引理了
我们要考虑的是不动点的数量，所以前 $gcd(i,n)$个点的连边方式确定下来，整个图的连边方式就确定了
我一开始想的是 $f[gcd(i,n)]$，但这是不对的
因为这些点也是可以看成分段的
最后一个点的下一个点可以看成是第一个点(即最后一段的下一段可以看成是第一段)，实际上这是一个环，所以这 $gcd(i,n)$个点的连接方式与 $F[gcd(i,n)]$是等价的，即合法性等效
所以答案是 $\sum^n_{i=1}F[gcd(i,n)]$
化简以后就是 $\sum^n_{i=1}[i|n]F(i)\varphi(\frac n i)$
万里长征就差最后一步了
除个n就可以了
啥？你说n,m可能不互质？
难道还要补个中国剩余定理……
有一个好用的式子
$a/b\mod p=(a\mod bp)/b$,其中 $b|a$
证明也十分简单
设 $a=bk$
$a/b\mod p=(bk)/b\mod p=k\mod p$
$(a\mod bp)/b=(bk\mod bp)/b=b(k\mod p)/b=k\mod p$
所以运算时模 $nm$，最后除以m就可以了
$nm≤10^{18}$所以做乘法的时候还要用快速乘……

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#define LL long long 
using namespace std;
int n,m;
int prime[32000];
bool vis[32000];
LL mo,ans;
LL mul(LL x,LL y)
{
    LL t=0;
    bool flag=0;
    if (y<0) y=-y,flag=1;
    for (;y;y>>=1,x=(x+x)%mo)
        if (y&1) t=(t+x)%mo;
    return flag?-t:t;
}
struct Matrix
{
    LL a[4][4];
    void clr(){memset(a,0,sizeof(a));}
}A,B;
Matrix operator *(Matrix A,Matrix B)
{
    Matrix C;
    C.clr();
    for (int i=1;i<=3;++i)
        for (int j=1;j<=3;++j)
            for (int k=1;k<=3;++k)
                C.a[i][j]=(C.a[i][j]+mul(A.a[i][k],B.a[k][j]))%mo;
    return C;
}
LL cal(int x)
{
    if (x==1) return 1;
    if (x==2) return 5;
    A.clr();B.clr();
    A.a[1][1]=5;A.a[1][2]=1;A.a[1][3]=2;
    B.a[1][1]=3;B.a[2][1]=-1;B.a[1][2]=1;
    B.a[3][1]=1;B.a[3][3]=1;
    x-=2;
    for (;x;x>>=1,B=B*B)
        if (x&1) A=A*B;
    return A.a[1][1];
}
int phi(int x)
{
    int ans=x;
    for (int i=1;prime[i]*prime[i]<=x;++i)
        if (x%prime[i]==0)
        {
            ans=ans/prime[i]*(prime[i]-1);
            while (x%prime[i]==0) x/=prime[i];
        }
    if (x>1) ans=ans/x*(x-1);
    return ans%mo;
}
main()
{
    int lim=31625;
    for (int i=2;i<=lim;++i)
    {
        if (!vis[i])
            prime[++prime[0]]=i;
        for (int j=1;j<=prime[0];++j)
        {
            if (i*prime[j]>lim) break;
            vis[i*prime[j]]=1;
            if (i%prime[j]==0) break;
        }
    }
    while (~scanf("%d%d",&n,&m))
    {
        mo=1LL*n*m;
        ans=0;
        for (int i=1;i*i<=n;++i)
        if (n%i==0)
        {
            ans=(ans+mul(phi(n/i),cal(i)))%mo;
            if (i*i!=n) ans=(ans+mul(phi(i),cal(n/i)))%mo;
        }
        printf("%d\n",(ans/n+m)%m);
    }
}