[BZOJ3456]城市规划：DP+NTT+多项式求逆

写在前面的话

昨天听吕老板讲课，数数题感觉十分的神仙。

于是，ErkkiErkko这个小蒟蒻也要去学数数题了。

分析

Miskcoo orz

带标号无向连通图计数。

(f(x))表示(x)个点的带标号无向连通图的个数。弱化限制条件，令(g(x))表示(x)个点的带标号无向图的个数（不要求连通）。

考虑每条边是否出现，显然有：

[g(x)=2^{inom{x}{2}} ]

考虑编号为(1)的结点所在连通块的大小，有：

[g(x)=sum_{i=1}^{x}inom{x-1}{i-1} imes f(i) imes g(x-i) ]

把组合数拆开，

[frac{g(x)}{(x-1)!}=sum_{i=1}^{x}frac{f(i)}{(i-1)!} imes frac{g(x-i)}{(x-i)!} ]

于是就可以用多项式求逆求出(f(n))了。

代码

#include <bits/stdc++.h>
#define rin(i,a,b) for(register int i=(a);i<=(b);++i)
#define irin(i,a,b) for(register int i=(a);i>=(b);--i)
#define trav(i,a) for(register int i=head[a];i;i=e[i].nxt)
typedef long long LL;
using std::cin;
using std::cout;
using std::endl;

inline int read(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}

const int MAXN=130005;
const int NTT=1048576;
const LL MOD=1004535809;
const LL G=3;
const LL INVG=334845270;

int N,n,m,len;
LL w[NTT+5],iw[NTT+5];
LL fac[MAXN],invf[MAXN];
LL g[MAXN];
LL rev[MAXN<<2],A[MAXN<<2],B[MAXN<<2];

inline LL qpow(LL x,LL y){
	LL ret=1,tt=x%MOD;
	while(y){
		if(y&1) ret=ret*tt%MOD;
		tt=tt*tt%MOD;
		y>>=1;
	}
	return ret;
}

void ntt(LL *c,int dft){
	rin(i,0,n-1)
		if(i<rev[i])
			std::swap(c[i],c[rev[i]]);
	for(register int mid=1;mid<n;mid<<=1){
		int r=(mid<<1),u=NTT/r;
		for(register int l=0;l<n;l+=r){
			int v=0;
			for(register int i=0;i<mid;++i,v+=u){
				LL x=c[l+i],y=c[l+mid+i]*(dft>0?w[v]:iw[v])%MOD;
				c[l+i]=x+y<MOD?x+y:x+y-MOD;
				c[l+mid+i]=x-y>=0?x-y:x-y+MOD;
			}
		}
	}
	if(dft<0){
		LL invn=qpow(n,MOD-2);
		rin(i,0,n-1) c[i]=c[i]*invn%MOD;
	}
}

void getinv(int mdx){
	if(mdx==1){
		A[0]=qpow(g[0],MOD-2);
		return;
	}
	getinv((mdx+1)>>1);
	m=(mdx-1)+((((mdx+1)>>1)-1)<<1),len=0;
	for(n=1;n<=m;n<<=1) ++len;
	rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
	rin(i,0,n-1) B[i]=i<mdx?g[i]:0;
	ntt(A,1);ntt(B,1);
	rin(i,0,n-1) A[i]=(2*A[i]-B[i]*A[i]%MOD*A[i]%MOD+MOD)%MOD;
	ntt(A,-1);
	rin(i,mdx,n-1) A[i]=0;
}

void init(){
	LL v=qpow(G,(MOD-1)/NTT),iv=qpow(INVG,(MOD-1)/NTT);
	w[0]=iw[0]=1;
	rin(i,1,NTT-1) w[i]=w[i-1]*v%MOD,iw[i]=iw[i-1]*iv%MOD;
	fac[0]=1;
	rin(i,1,N) fac[i]=fac[i-1]*i%MOD;
	invf[N]=qpow(fac[N],MOD-2);
	irin(i,N-1,0) invf[i]=invf[i+1]*(i+1)%MOD;
}

int main(){
	N=read();
	if(N==0){
		printf("1
");
		return 0;
	}
	init();
	g[0]=1;
	rin(i,1,N) g[i]=qpow(2,1ll*i*(i-1)/2)*invf[i]%MOD;
	getinv(N+1);
	m=(N<<1),len=0;
	for(n=1;n<=m;n<<=1) ++len;
	rin(i,0,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
	B[0]=0;
	rin(i,1,N) B[i]=qpow(2,1ll*i*(i-1)/2)*invf[i-1]%MOD;
	ntt(A,1);ntt(B,1);
	rin(i,0,n-1) A[i]=A[i]*B[i]%MOD;
	ntt(A,-1);
	printf("%lld
",A[N]*fac[N-1]%MOD);
}