HDU 6159 Graph Of Zhuper（dp+CDQ分治+FFT+生成函数+莫比乌斯反演）

Description

$G$为一个$n$个点带标号无向图，无重边无自环，假设其有$k$个连通分支，第$i$个连通分支点数为$size_i$，定义$G$的权值$Zhu(G)=gcd(size_1,size_2,...,size_k)$，求$\sum Zhu(G)$

Input

第一行一整数$T$表示用例组数，每组用例输入两个整数$n,P$分别表示图的点数和模数$(1\le T\le 10,1\le n\le 30000,10^8\le P\le 10^9+9)$

Output

对于每组用例，输出答案，结果模$P$

Sample Input

3
1 1000000007
2 1000000007
3 1000000007

Sample Output

Case #1: 1
Case #2: 3
Case #3: 16

Solution

令$dp[n]$为$n$个点带标号连通无向图个数，$a[n]$为$n$个点带标号无向图的个数，则$a[n]=2^{\frac{n(n-1)}{2}}$，为求$dp[n]$，枚举第一个点所在连通块点数$k$，则有$dp[n]=a[n]-\sum\limits_{k=1}^{n-1}C_{n-1}^{k-1}\cdot dp[k]\cdot a[n-k]$，即$\frac{dp[n]}{(n-1)!}=n\cdot \frac{a[n]}{n!}-\sum\limits_{k=1}^{n-1}\frac{dp[k]}{(k-1)!}\frac{a[n-k]}{(n-k)!}$，令$A[n]=\frac{dp[n]}{(n-1)!},B[n]=\frac{a[n]}{n!}$，下面考虑用$CDQ$分治求$A[n]$，假设已经求出$A[l],A[l+1],...,A[mid]$，考虑求出这些值对$A[mid+1],A[mid+2],...,A[r]$的影响，对于$t\in[mid+1,r]$，设$A[l],A[l+1],...,A[mid]$对其影响为$-C(t)$，则有$C(t)=\sum\limits_{k=l}^{mid}A[k]\cdot B[n-k]$，$FFT$即求出影响序列$C$，进而得到$dp$序列

令$h(d)$为$n$个点带标号无向图，每个连通块点数被$d$整除的个数，则$d|n$，先不考虑点的标号，最后乘上$n!$即可，则一个点数被$d$整除的连通块的生成函数为$f(x)=\sum\limits_{i=1}^{n/d}\frac{dp[id]}{(id)!}x^i$，多个连通块的生成函数$g(x)=\sum\limits_{i=0}^{\infty}\frac{f^i(x)}{i!}=e^{f(x)}$，做多项式$exp$即可得到$g(x)$，进而$h(d)=n!\cdot [x^{n/d}]g(x)$

令$H(d)$为$n$个点带标号无向图，连通块点数的最大公约数为$d$的个数，则有$h(d)=\sum\limits_{d|m}H(m)$，由莫比乌斯反演，$H(d)=\sum\limits_{d|m}\mu(\frac{m}{d})h(m)$， 答案即为$ans=\sum\limits_{d|n}d\cdot H(d)=\sum\limits_{d|n}d\cdot \sum\limits_{d|m}\mu(\frac{m}{d})\cdot h(m)=\sum\limits_{m|n}h(m)\cdot \sum\limits_{d|m}d\cdot \mu(\frac{m}{d})=\sum\limits_{m|n}h(m)\cdot \varphi(m)$

Code

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;
#define maxn 30005
#define maxfft 131072+5
int mod;
const double pi=acos(-1.0);
struct cp 
{
    double a,b;
    cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
    cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
    cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
    cp operator *(const double &o)const {return (cp){a*o,b*o};}
    cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
    int j=0;
    while((1<<j)<len)j++;
    j--;
    for(int i=0;i<len;i++)
        pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
    for(int i=0;i<len;i++)
        if(i<pos[i])swap(x[i],x[pos[i]]);
    w[0]=(cp){1,0};
    for(unsigned i=2;i<=len;i<<=1)
    {
        cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
        for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
        for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
        for(int j=0;j<len;j+=i)
        {
            cp *a=x+j,*b=a+(i>>1);
            for(int l=0;l<i>>1;l++)
            {
                cp o=b[l]*w[l];
                b[l]=a[l]-o;
                a[l]=a[l]+o;
            }
        }
    }
    if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
int temp[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
    if(n<=100&&m<=100||min(n,m)<=5)
    {
        for(int i=0;i<n+m-1;i++)temp[i]=0;
        for(int i=0;i<n;i++)
            for(int j=0;j<m;j++)
            {
                temp[i+j]+=(ll)a[i]*b[j]%mod;
                if(temp[i+j]>=mod)temp[i+j]-=mod;
            }
        for(int i=0;i<n+m-1;i++)c[i]=temp[i];
        return ;
    }
    int len=1;
    while(len<n+m)len<<=1;
    fft_init(len);
    for(int i=0;i<len;i++)
    {
        int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
        x[i]=(cp){(aa>>15),(aa&32767)},y[i]=(cp){(bb>>15),(bb&32767)};
    }
    fft(x,len,1),fft(y,len,1);
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod;
        ta=(ta<<15)%mod;
        c[i]=ta;
    }
    for(int i=0;i<len;i++)
    {
        int j=len-1&len-i;
        z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
    }
    fft(z,len,-1);
    for(int i=0;i<n+m-1;i++)
    {
        ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
        ta=(ta+(tb<<30))%mod;
        c[i]=(c[i]+ta)%mod;
    }
}
int mod_pow(int a,ll b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=(ll)ans*a%mod;
        a=(ll)a*a%mod;
        b>>=1;
    }
    return ans;
}
int T,n,a[maxn],b[maxfft],dp[maxn],fact[maxn],ifact[maxn],inv[maxn],phi[maxn];
void init(int n)
{
    fact[0]=1;
    for(int i=1;i<=n;i++)fact[i]=(ll)i*fact[i-1]%mod;
    inv[1]=1;
    for(int i=2;i<=n;i++)inv[i]=mod-(ll)(mod/i)*inv[mod%i]%mod;
    ifact[0]=1;
    for(int i=1;i<=n;i++)ifact[i]=(ll)inv[i]*ifact[i-1]%mod;
    a[0]=0;
    for(int i=1;i<=n;i++)a[i]=(ll)mod_pow(2,(ll)i*(i-1)/2)*ifact[i]%mod;
    for(int i=1;i<=n;i++)phi[i]=i;
    for(int i=2;i<=n;i++)
        if(phi[i]==i)
            for(int j=i;j<=n;j+=i)
                phi[j]=phi[j]/i*(i-1);
}
void deal(int l,int r)
{
    if(l==r)
    {
        dp[l]+=(ll)l*a[l]%mod;
        if(dp[l]>=mod)dp[l]-=mod;
        return ; 
    }
    int mid=(l+r)>>1;
    deal(l,mid);
    FFT(dp+l,a+1,mid-l+1,r-l,b);
    for(int i=mid+1;i<=r;i++)
    {
        dp[i]-=b[i-l-1];
        if(dp[i]<0)dp[i]+=mod;
    }
    deal(mid+1,r);
}
int temp1[maxfft],temp2[maxfft],temp3[maxfft];
void Poly_Inv(int *poly,int n,int *ans)
{
    ans[0]=1;
    for(int i=2;i<=n;i<<=1)
    {
        FFT(poly,ans,i,i/2,temp1);
        FFT(ans,temp1+i/2,i/2,i/2,temp1);
        for(int j=0;j<i/2;j++)ans[j+i/2]=mod-temp1[j];
    }
}
void Poly_Log(int *poly,int n,int *ans)
{
    Poly_Inv(poly,n,temp2);
    for(int i=0;i<n-1;i++)ans[i]=(ll)poly[i+1]*(i+1)%mod;
    FFT(ans,temp2,n-1,n,ans);
    for(int i=n-1;i>0;i--)ans[i]=(ll)ans[i-1]*inv[i]%mod;
    ans[0]=0;
}
void Poly_Exp(int *poly,int n,int *ans)
{
    if(n==1)
    {
        ans[0]=1;
        return ;
    }
    Poly_Exp(poly,n/2,ans);
    Poly_Log(ans,n,temp3);
    for(int i=0;i<n;i++)
    {
        temp3[i]=poly[i]-temp3[i];
        if(temp3[i]<0)temp3[i]+=mod;    
    }
    temp3[0]++;
    if(temp3[0]==mod)temp3[0]=0;
    FFT(ans,temp3,n,n,ans);
    for(int i=n;i<2*n;i++)ans[i]=0;
}
int f[maxfft],g[maxfft];
int Solve(int n,int d)
{
    if(d==1)return mod_pow(2,(ll)n*(n-1)/2);
    int m=n/d+1;
    f[0]=0;
    for(int i=1;i<m;i++)f[i]=(ll)dp[i*d]*ifact[i*d]%mod;
    int len=1;
    while(len<m)len<<=2;
    for(int i=m;i<len;i++)f[i]=0;
    Poly_Exp(f,len,g);
    return (ll)g[n/d]*fact[n]%mod;
}
int main()
{
    int Case=1;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&n,&mod);
        init(n);
        memset(dp,0,sizeof(dp));
        deal(1,n);
        for(int i=1;i<=n;i++)dp[i]=(ll)fact[i-1]*dp[i]%mod;
        int ans=0;
        for(int i=1;i*i<=n;i++)
            if(n%i==0)
            {
                ans+=(ll)Solve(n,i)*phi[i]%mod;
                if(ans>=mod)ans-=mod;
                if(i*i!=n)
                {
                    ans+=(ll)Solve(n,n/i)*phi[n/i]%mod;
                    if(ans>=mod)ans-=mod;
                }
            }
        printf("Case #%d: %d\n",Case++,ans);
    }
    return 0;
}