bzoj3456/jzoj3303：城市规划（画柿子+多项式逆元+NTT）

我noip初赛最后半个小时就在数4个点的连通图有几个，结果还数错了…然后我就很膜拜这题，终于被我在jzoj找到了。

权限题传送门

题意：n个不同的点的构成的连通图有多少种，模一个费马素数。n≤2e5

设f[i]为i个点的连通图方案数，思路为所有方案减去不合法的。

所有方案数为$2^{i*(i-1)/2}$。
对于不连通的，枚举1号点所在连通块大小为j。选这剩下j-1个点的方案数为$C_{i-1}^{j-1}$。连通块内方案数为f[j]，连通块外方案数为$2^{(i-j)*(i-j-1)/2}$。

整理一下，有

$$f[i]=2^{i*(i-1)/2}-\sum_{j=1}^{i-1}f[j]*C_{i-1}^{j-1}*2^{(i-j)*(i-j-1)/2}$$
两边同除(i-1)!，有 $$\frac{f[i]}{(i-1)!}=\frac{2^{i*(i-1)/2}}{(i-1)!}-\sum_{j=1}^{i-1}\frac{f[j]*2^{(i-j)*(i-j-1)/2}}{(j-1)!*(i-j)!}$$

想到若1号点所在的连通块大小为i的含义，尝试将和式的j取i，发现它正好等于左边，有

$$\sum_{j=1}^i\frac{f[j]*2^{(i-j)*(i-j-1)/2}}{(j-1)!*(i-j)!}=\frac{2^{i*(i-1)/2}}{(i-1)!}$$
继续 $$\sum_{j=1}^i\frac{f[j]}{(j-1)!}*\frac{2^{(i-j)*(i-j-1)/2}}{(i-j)!}=\frac{2^{i*(i-1)/2}}{(i-1)!}$$

左边已经是个卷积了。我们要求的f在左边，求个多项式逆元就可以了。

老是T的NTT

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>

#define mmst(a, b) memset(a,b,sizeof(a))
#define mmcp(a, b) memcpy(a,b,sizeof(a))

using namespace std;

typedef long long LL;

const int N=400400;
const LL p=1004535809,g=3;

LL cheng(LL a,LL b)
{
    LL res=1;
    for(;b;b>>=1,a=a*a%p)
    if(b&1)
    res=res*a%p;
    return res;
}

int n,rev[N];

void init(int lim)
{
    n=1;
    int k=-1;
    while(n<lim)
    n<<=1,k++;
    for(int i=0;i<n;i++)
    rev[i]=(rev[i>>1]>>1) | ((i&1)<<k);
}

void ntt(LL *a,bool ops)
{
    for(int i=0;i<n;i++)
    if(i<rev[i])
    swap(a[i],a[rev[i]]);

    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        LL wn;
        if(ops)
        wn=cheng(g,(p-1)/l);
        else
        wn=cheng(g,p-1-(p-1)/l);
        for(int i=0;i<n;i+=l)
        {
            LL w=1;
            for(int k=0;k<m;k++)
            {
                LL t=a[i+k+m]*w%p;
                a[i+k+m]=(a[i+k]-t+p)%p;
                a[i+k]=(a[i+k]+t)%p;
                w=w*wn%p;
            }
        }
    }
    if(!ops)
    {
        LL Inv=cheng(n,p-2);
        for(int i=0;i<n;i++)
        a[i]=a[i]*Inv%p;
    }
}

int nn;
LL jc[N],Ijc[N],B[N],C[N],D[N],E[N];

void ny(int x)
{
    if(x==1)
    return;
    ny(x/2);

    for(int i=0;i<x;i++)
    E[i]=B[i];
    init(x*2);
    ntt(E,1);
    ntt(D,1);
    for(int i=0;i<n;i++)
    D[i]=D[i]*(2-E[i]*D[i]%p+p)%p;
    ntt(D,0);

    for(int i=0;i<n;i++)
    E[i]=0;
    for(int i=x;i<n;i++)
    D[i]=0;
}

int main()
{
    cin>>nn;
    jc[0]=1;
    for(int i=1;i<=nn;i++)
    jc[i]=jc[i-1]*i%p;

    for(int i=0;i<=nn;i++)
    Ijc[i]=cheng(jc[i],p-2);

    for(LL i=1;i<=nn;i++)
    C[i]=cheng(2,i*(i-1)/2)*Ijc[i-1]%p;

    for(LL i=0;i<=nn;i++)
    B[i]=cheng(2,i*(i-1)/2)*Ijc[i]%p;
    B[0]=1;

    D[0]=1;

    init(nn+2);
    ny(n);

    init(nn*2+2);
    ntt(D,1);
    ntt(C,1);
    for(int i=0;i<n;i++)
    C[i]=C[i]*D[i]%p;
    ntt(C,0);

    cout<<C[nn]*jc[nn-1]%p<<endl;

    return 0;
}