【数论-什么都有的数论】BZOJ3113（HDU2481）Toy

【前言】
这是一道什么都有的数论题，我想的时候一脸懵逼，根本不会做。
在使用“特殊方法”后搞了出来，但是码的过程也一波三折。
一道神题！

【题目】
BZOJ地址
HDU地址
外面有一圈$n$个结点，中心有一个结点与$n$个结点都相连，总共就是$2\times n$条边，删除$n$条边，使$n+1$个点连通，旋转相同视为等价，问有多少种情况。

【分析】
我们接下来将看到：素数筛选，求解欧拉函数，BurnSide引理，二分模拟乘法（快速乘），递推的构造，矩阵快速幂，置换群，枚举……

【解题思路】
参考博文

根据分析发现旋转似乎并不重要，我们算出所有结果后除以$n$就行了。
我们考虑任取两个节点$a,b$讨论。

那么总数便是$a,b$断开的种数与$a,b$连在一起的种数的和。
$f(n)$表示外圈有$n$个结点时，而$a,b$是断开的种数。
$g(n)$表示外圈有$n$个结点时，而$a,b$是连在一起的种数。
如果$a,b$之间是断开的，此时若与$a$直接相连的为$k$个(加上$a$自己)，那么显然这$k$个要与其它的保持连通的，与中心必须有一条边，如果有多条边就形成环了，显然不满足生成树。另外$n-k$个点就为$f[n-k]$种，我们可以枚举$k$，则$f[n]=\sum (i\times f[n-i]) (n-1\geq i \geq 0)$

如果$a,b$是连在一起的，此时若与$a,b$相连的为$k$个(包括$a,b$)，那么$a,b$是相邻的在这$k$个位置选择就有$k-1$种，而这$k$个与中心相连的选择有$k$种，剩下的与这部分是分开的，则为$f[n-k]$，所以可以枚举$k$，最终结果$g[n]=\sum (i\times (i-1)\times f[n-i]) (n-1\geq i\geq 2)$

则最终的种数便是$T[n]=f[n]+g[n]$。

接下来就是式子的化简，这部分我倒是自己尝试了一下，不算很难。

$$\begin{eqnarray} f[n]&=&\sum (i\times f[n-i]) (n-1\geq i \geq 0) \\ &=& f[n-1]+2\times f[n-2]+3\times f[n-3]\dots (n-1)\times f[1]+n\times f[0]\\ &=& f[n-1]+f[n-2]\dots f[1]+f[0]+(f[n-2]+2\times f[n-3]\dots +(n-1)\times f[0]) \end{eqnarray}$$
令 $s[n]$为 $f[i]$的前 $n$项和，则有：
$$\begin{eqnarray} f[n]&=&s[n-1]+f[n-2]+2\times f[n-3]\dots (n-2)f[1] \\ &=& s[n-1]+\sum(i*f((n-1)-i)) (n-2\geq i\geq 0)\\ &=& s[n-1]+f[n-1]\\ &=& s[n-2]+f[n-1]+f[n-1]\\ &=& 3\times f[n-1]+f[n-2]\cdots \cdots(\because s[n-1]=f[n]-f[n-1]) \end{eqnarray}$$

差不多的我们继续推$g[n]$

$$\begin{eqnarray} g[n]&=&\sum (i\times (i-1)\times f[n-i]) (n-1\geq i\geq 2) \\ &=&1\times2\times f[n-2]+2\times3\times f[n-3]+3\times4\times f[n-4]\dots (n-1)\times (n-2)\times f[1]\\ \end{eqnarray}$$

那么我们有：

$$\begin{eqnarray} g[n-1]&=&1\times 2\times f[n-3]+2\times 3\times f[n-4]+\dots +(n-2)\times (n-3)\times f[1] \\ \end{eqnarray}$$

这样：

$$\begin{eqnarray} g[n]-g[n-1]&=&2\times f[n-2]+4\times f[n-3]+\dots +(2\times (n-2))*f[1]\\ &=&2\times (f[n-2]+2*f[n-3]+\dots +f[1])\\ &=&2\times f[n-1] \end{eqnarray}$$

因此：

$$\begin{eqnarray} g[n]&=&2\times (f[1]+f[2]+f[3]+\dots +f[n-1])\\ &=&2\times (s[n-1]-f[0])\\ &=&2\times(f[n]-f[n-1]-1) \end{eqnarray}$$

其中$f[0]=1$，但$s[n]$是不包含$f[0]$的
处理到这里，我们可以用矩阵快速幂求出$f[n]$，顺便也可以将$g[n]$搞出来。
所以$T[n]$就解决咯。

关于矩阵的构造，我们令：

$$A= \left( \begin{array}{cc} 3 & 1 \\ -1 & 0 \end{array} \right)$$
那么可以得到：
$$\left( \begin{array}{cc} f[n] & f[n-1] \end{array} \right)= \left( \begin{array}{cc} f[1] & f[0] \end{array} \right)\times A^{n-1}$$

接下来还要上Burnside定理+欧拉函数优化。

我们要考虑的是不动点的数量，所以前$gcd(i,n)$个点的连边方式确定下来，整个图的连边方式就确定了
最后一个点的下一个点可以看成是第一个点(即最后一段的下一段可以看成是第一段)，实际上这是一个环，所以这$gcd(i,n)$个点的连接方式与$T[gcd(i,n)]$是等价的，即合法性等效
所以答案是$\sum^n_{i=1}T[gcd(i,n)]$
化简以后就是$\sum^n_{i=1}(T[i]\times \phi (\frac n i))(i|n)$

大概应该是这样的吧，可能最后一步有点不对qwq。
然后因为$n,m$不一定互质，我们也许还要用CRT。
但其实我们可以用这个式子：

$$\frac a b \mod p=\frac{(a \mod bp)} b,其中b|a$$
证明也十分简单
设 $a=bk$
$$\frac a {b \mod p}=\frac {bk} {b\mod p}=k \mod p$$
$$\frac {a \mod bp} b=\frac {bk \mod bp} b=\frac {b(k \mod p)} b= k\mod p$$
所以运算时模 $nm$，最后除以 $m$就可以了
$nm\leq 10^{18}$所以做乘法的时候还要用快速乘

以上，就是这道神题了。

接下来是作弊时间
暴力打出前六个数字：1，3，6，13，25，58
丢进OEIS
得到规律：1/n*Sum_{d divides n} phi(n/d)*A004146(d)
A004146的公式是这样： a(n+1) = 3*a(n) - a(n-1) + 2

以上就是一种珂学的方法了。

【参考代码】

#include<bits/stdc++.h>
using namespace std;

typedef long long LL;
const int N=4e4+4;
const int M=1e9+10;
int n,cnt;
LL mod;
bool flag[N];
int pri[N];

struct Matrix
{
    LL m[2][2];

    Matrix(){
        memset(m,0,sizeof(m));
    }
};
Matrix init;

LL mul(LL a,LL b)
{
    a%=mod;b%=mod;
    if(a<0) a+=mod;
    if(b<0) b+=mod;
    LL ret=0;
    while(b)
    {
        if(b&1)
        {
            ret+=a;
            if(ret>=mod)
                ret-=mod;
        }
        a<<=1ll;
        if(a>=mod)
            a-=mod;
        b>>=1ll;
    }
    return ret;
}

Matrix operator * (Matrix a,Matrix b)
{
    Matrix ret;
    for(int i=0;i<2;++i)
        for(int j=0;j<2;++j)
            for(int k=0;k<2;++k)
                (ret.m[i][j]+=mul(a.m[i][k],b.m[k][j]))%=mod;
    return ret;
}

Matrix operator ^ (Matrix a,int b)
{
    Matrix ret;
    for(int i=0;i<2;++i)
        for(int j=0;j<2;++j)
            ret.m[i][j]=(i==j);
    while(b)
    {
        if(b&1)
            ret=ret*a;
        a=a*a;
        b>>=1;
    }
    return ret;
}

void Prime()
{
    int t=sqrt(M);
    for(int i=0;i<=t;++i)
        flag[i]=false;
    flag[1]=true;

    for(int i=2;i<=t;++i)
    {
        if(!flag[i])
            pri[cnt++]=i;
        for(int j=0;j<cnt && i*pri[j]<=t;++j)
        {
            flag[i*pri[j]]=true;
            if(!(i%pri[j]))
                break;
        }
    }
}

int Eular(int t)
{
    int ret=1;
    for(int i=0;i<cnt && pri[i]*pri[i]<=t;++i)
    {
        if(!(t%pri[i]))
        {
            t/=pri[i];ret*=pri[i]-1;
            while(!(t%pri[i]))
            {
                t/=pri[i];
                ret=(ret*pri[i])%mod;
            }
        }
    }
    if(t>1)
        ret*=t-1;
    return ret%mod;
}

LL Get_T(int k)
{
    if(k==1)    
        return 1;
    else
    if(k==2)
        return 5;
    Matrix tmp=init^(k-2);
    LL f=3ll*tmp.m[0][0]+1ll*tmp.m[1][0]%mod;
    LL g=2ll*(f-(3ll*tmp.m[0][1]+tmp.m[1][1])-1)%mod;
    return ((g+f)%mod+mod)%mod;
}

LL Polya()
{
    LL sum=0;
    int i;
    for(i=1;i*i<n;++i)
        if(!(n%i))
        {
            sum=(sum+mul(Eular(i),Get_T(n/i)))%mod;
            sum=(sum+mul(Eular(n/i),Get_T(i)))%mod;
        }
    if(i*i==n)
        sum=(sum+mul(Get_T(i),Eular(i)))%mod;
    return sum/n;
}

int main()
{
//  freopen("BZOJ3113.in","r",stdin);
//  freopen("BZOJ3113.out","w",stdout);

    Prime();
    init.m[0][0]=3;init.m[0][1]=1;init.m[1][0]=-1;init.m[1][1]=0;
    while(scanf("%d%lld",&n,&mod)!=EOF)
    {
        mod=1ll*n*mod;
        printf("%lld\n",Polya()%(mod/n));
    }

    return 0;
}

【总结】
数论真是太可怕了！我们还是OEIS吧！