[BZOJ3456] 城市规划 - 快速傅里叶变换 - 快速数论变换 - 卷积 - 多项式求逆

-太穷买不起权限号只能嘴巴+其他网站测评AC

----------以下部分为纸张蒟蒻zzt想到的部分，神犇可忽略不计----------

    啥子？要求代标号无向连通图？这不是思博题？不是一眼就转化为图总数-不连通的数目？蛤？不连通数目是啥？QAQ窝觉得是∑大小为i的连通图数目*大小为n-i的非连通图数目…（旁边：这不是一脸错误？泥算算小数据对吗昂）太弱被狂D。

----------以上部分就是纸张蒟蒻zzt想到的显而易见的性质，神犇可忽略不计----------

题解：

设g(i)=大小为i的图总数=2^(i*(i-1)/2)，f(i)=大小为i的连通图数目

枚举第一个点的连通块大小，则g(n)=∑C(n-1,j-1)f(j)g(n-j),j=1...n

显然我们要求的是f，那么我们怎么求f(n)呢

将g(i)的表达式带入原式，得2^(n*(n-1)/2)=∑f(j)*(n-1)!/(j-1)!/(n-j)!*2^((n-j)*(n-j-1)/2)

那么将两边同处以(n-1)!，发现式子对应生成函数的卷积，令A=∑2^(i*(i-1)/2)/(i-1)!*x^i，B=∑2^(i*(i-1)/2)/i!*x^i，C=∑f(i)*x^i

得A=BC 即C=A*(B^-1)

那么我们只要求出B在mod x^(n+1)意义下的逆元再用NTT与A相乘即可。

求逆元可以分块FFT，先求Bmod x^floor(n/2)意义下的逆元为S'，则Bmod x^(n+1)意义下的逆元S满足：

BS≡1 BS'≡1(mod x^floor(n/2)) 即B(S-S')≡0，显然B!=0(i=0...∞) 所以S-S'≡0(mod x^(floor(n/2))

即S^2-2SS'+S'^2≡0(mod x^(n+1))

同乘B降次，得S=2S'-BS'^2，后面的部分NTT即可

#include"bits/stdc++.h"

using namespace std;
typedef long long ll;

const ll P=(479<<21)+1,g=3;
const int N=262150,LG=20;
ll w[LG];

ll ksm(ll a,ll t,ll mod){
	ll re=1;
	while(t){
		if(t%2)re=re*a%mod;
		a=a*a%mod;t>>=1;
	}
	return re;
}

void build(int len){
	for(int i=1;(1<<i)<=len;i++)
		w[i]=ksm(g,(P-1)>>i,P);
}

template <class __t>
void Rader(__t a[],int len){
	int i,j=len>>1,k;
	for(i=1;i<len-1;i++){
		if(i<j)swap(a[i],a[j]);
		for(k=len>>1;j>=k;j-=k,k>>=1);
		if(j<k)j+=k;
	}
}

void NTT(ll a[],int len,int DFT){
	int i,j,k,l,t;Rader(a,len);
	for(t=1;(1<<t)<=len;t++){
		l=1<<t;i=l>>1;
		for(j=0;j<len;j+=l){
			ll wn=1;
			for(k=j;k<j+i;k++){
				ll u=a[k],v=a[k+i]*wn%P;
				a[k]=(u+v)%P;a[k+i]=(u-v+P)%P;
				wn=wn*w[t]%P;
			}
		}
	}
	if(DFT==-1){
		for(i=1;i<len>>1;i++)
			swap(a[i],a[len-i]);
		ll mul=ksm(len,P-2,P);
		for(i=0;i<len;i++)
			a[i]=a[i]*mul%P;
	}
}

void Convol(ll a[],ll b[],ll in[],int len){
	NTT(a,len,1);NTT(b,len,1);
	for(int i=0;i<len;i++)
		in[i]=a[i]*b[i]%P;
	NTT(in,len,-1);
}

void SQR(ll a[],int len){
	NTT(a,len,1);
	for(int i=0;i<len;i++)
		a[i]=a[i]*a[i]%P;
	NTT(a,len,-1);
}

ll s[N],x[N],y[N];
void inverse(ll a[],int deg){
	if(deg==1){s[0]=ksm(a[0],P-2,P);return;}
	int mid=(deg+1)>>1,i,L;inverse(a,mid);
	memcpy(x,s,sizeof(s));memset(y,0,sizeof(y));
	for(i=0;i<deg;i++)y[i]=a[i];
	for(L=1;L<=mid*2;L<<=1);SQR(x,L);
	for(;L<=deg*2;L<<=1);Convol(x,y,x,L);
	for(i=0;i<deg;i++)s[i]=(2*s[i]-x[i]+P)%P;
}

ll F[N],G[N],ans[N];
ll fac[N>>1],invfac[N>>1],inv[N>>1];
int pri[N>>1],flag[N>>1],cnt,n,l;

void prepare(){
	scanf("%d",&n);int i,j;
	for(l=1;l<=(n+1)<<1;l<<=1);build(l);
	for(fac[1]=invfac[1]=inv[1]=1,i=2;i<=n;i++){
		fac[i]=fac[i-1]*i%P;
		if(flag[i]==0){
			pri[++cnt]=i;
			inv[i]=ksm(i,P-2,P);
		}
		for(j=1;pri[j]*i<=n;j++){
			inv[pri[j]*i]=inv[pri[j]]*inv[i]%P;
			if(i%pri[j]==0) break;
		}
		invfac[i]=invfac[i-1]*inv[i]%P;
	}
	F[1]=G[1]=G[0]=1;
	for(i=2;i<=n;i++){
		ll kuai=ksm(2,i-1,P);
		G[i]=G[i-1]*kuai%P*inv[i]%P;
		F[i]=F[i-1]*kuai%P*inv[i-1]%P;
	}
}

void work(){
	inverse(G,n+1);
	Convol(F,s,ans,l);
	printf("%lld\n",ans[n]*fac[n-1]%P);
}

int main(){
	prepare();work();
	return 0;
}