bzoj3456 城市规划【生成函数+多项式求逆元】

解题思路

我们设$f(n)$表示$n$个点组成有编号无向连通图的方案数，但不是很好直接求。

那再设$g(n)$表示$n$的点组成有编号无向图的方案数，由于无向图是由多个连通子图构成，那么为了避免重复计数，我们枚举1号点所在连通块大小，其他边随便连，即可得到dp方程：

$g(n)=\sum\limits_{i=1}^nC_{n-1}^{i-1}f(i)*g(n-i)$

而容易得到$g(n)=2^{C_n^2}$，代入上式得：

$2^{C_n^2}=\sum\limits_{i=1}^nC_{n-1}^{i-1}f(i)*2^{C_{n-i}^2}$

两边同除以$(n-1)!$，则有：

$\frac{2^{C_n^2}}{(n-1)!}=\sum\limits_{i=1}^n\frac{f(i)}{(i-1)!}*\frac{2^{C_{n-i}^2}}{(n-i)!}$

发现这是一个卷积形式，考虑生成函数：

$F(x)=\sum\limits_{n=0}^\infty\frac{f(n)}{(n-1)!}$

$C(x)=\sum\limits_{n=0}^\infty\frac{2^{C_n^2}}{(n-1)!}$

$G(x)=\sum\limits_{n=0}^\infty\frac{2^{C_n^2}}{n!}$

则有：

$C(x)=F(x)G(x)$

$F(x)=C(x)G^{-1}(x)$

所以我们只需求出$G(x)$在$mod(x^{n+1})$意义下的逆元即可。

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

int getint()
{
    int i=0,f=1;char c;
    for(c=getchar();c!='-'&&(c<'0'||c>'9');c=getchar());
    if(c=='-')f=-1,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
    return i*f;
}

const int N=500005,p=1004535809,g=3;
int n,pos[N],a[N],b[N],c[N],w[N],w_inv[N],fac[N],fac_inv[N],bin[N];

int Pow(int x,ll y)
{
    int res=1;
    for(;y;y>>=1,x=(ll)x*x%p)
        if(y&1)res=(ll)res*x%p;
    return res;
}

void NTT(int f[],int len,int on)
{
    for(int i=1;i<len;i++)
        pos[i]=(i&1)?pos[i>>1]>>1|(len>>1):pos[i>>1]>>1;
    for(int i=1;i<len;i++)
        if(i<pos[i])swap(f[i],f[pos[i]]);
    for(int i=1,num=1;i<len;i<<=1,num++)
    {
        int wn=(on==1?w[num]:w_inv[num]);
        for(int j=0;j<len;j+=(i<<1))
        {
            int wi=1;
            for(int k=j;k<j+i;k++)
            {
                int u=f[k],v=(ll)wi*f[k+i]%p;
                f[k]=(u+v)%p,f[k+i]=(u-v+p)%p;
                wi=(ll)wi*wn%p;
            }
        }
    }
    if(on==-1)
        for(int i=0;i<len;i++)
            f[i]=(ll)f[i]*Pow(len,p-2)%p;
}

void Poly_inv(int a[],int b[],int deg)
{
    if(deg==1)
    {
        b[0]=Pow(a[0],p-2);
        return;
    }
    Poly_inv(a,b,(deg+1)>>1);
    static int tmp[N];
    int len=1;
    while(len<(deg<<1))len<<=1;
    for(int i=0;i<deg;i++)tmp[i]=a[i];
    for(int i=deg;i<len;i++)tmp[i]=0;
    for(int i=(deg+1)>>1;i<len;i++)b[i]=0;
    NTT(tmp,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++)b[i]=(ll)b[i]*((2-(ll)tmp[i]*b[i]%p+p)%p)%p;
    NTT(b,len,-1);
}

void multi(int f1[],int f2[])
{
    int len=1;
    while(len<(n<<1))len<<=1;
    NTT(f1,len,1),NTT(f2,len,1);
    for(int i=0;i<len;i++)f1[i]=(ll)f1[i]*f2[i]%p;
    NTT(f1,len,-1);
}

int main()
{
    //freopen("lx.in","r",stdin);
    //freopen("lx.out","w",stdout);
    n=getint();
    int len=1,num=0;
    while(len<((n+1)<<1))
    {
        len<<=1;
        w[++num]=Pow(g,(p-1)/len);
        w_inv[num]=Pow(w[num],p-2);
    }
    fac[0]=fac_inv[0]=1;
    for(int i=1;i<=n;i++)fac[i]=(ll)fac[i-1]*i%p;
    fac_inv[n]=Pow(fac[n],p-2);
    for(int i=n-1;i>=1;i--)fac_inv[i]=(ll)fac_inv[i+1]*(i+1)%p;
    bin[0]=bin[1]=1;
    for(int i=2;i<=n;i++)bin[i]=Pow(2,(ll)i*(i-1)/2);
    for(int i=0;i<=n;i++)a[i]=(ll)bin[i]*fac_inv[i]%p;
    for(int i=1;i<=n;i++)b[i]=(ll)bin[i]*fac_inv[i-1]%p;
    Poly_inv(a,c,n+1);
    for(int i=n+1;i<len;i++)c[i]=0;
    multi(c,b);
    int ans=(ll)c[n]*fac[n-1]%p;
    printf("%d\n",ans);
    return 0;
}