【分治NTT/多项式求逆】JZOJ3303. 城市规划

城市规划

• Description
• Solution
• 分治NTT
• 多项式求逆
• 分治NTT的代码
• 多项式求逆的代码
• 注意一个细节

Description

求出n 个点的简单(无重边无自环)无向连通图数目.
n <= 130000

Solution

• 设 f[i] 表示大小为i的答案，g[i]为2^C(i,2)表示大小为i的任意无向图个数。
• 考虑运用容斥。f[i] = g[i]-sigma( f[j] * g[i-j] * C(i-1,j-1) )
• 用所有方案减去不连通的方案。枚举1所在的连通块（保证不重复不遗漏），再考虑哪些点在这个联通块中，最后剩下的点随便连（反正不跟1所在联通块连），那么这个方案就是不连通的方案数。

分治NTT

• 将组合数拆开即可得到
  f[i] = g[i] - sigma( f[j] * g[i-j] * (i-1)! / (j-1)! / (i-j)!)

• 我们不难发现sigma里面有的只跟i-j有关，有的只跟j有关，所以我们将这些项整理一下，将无关的(i-1)!提出来，得到：
  f[i] = g[i] - (i-1)! * sigma( f[j] / (j-1)! * g[i-j] / (i-j) !)

• 设a[i]=f[i]/(i-1)! , b[i]=g[i]/i!

• sigma中的东西可以表示成a与b的卷积。

• 卷积的一般形式：F[ n ] = sigma (a[ i ] * b [ j ] | ( i + j = n ) )

• 现在我们已经将它转化成NTT的形式了，但是每一个i都是由前面全部的j所贡献得到的，我们不可能枚举每一个点再枚举它对后面的贡献，我们考虑分治。

• 当前区间[l,r]，先做[l,mid]，再考虑左边对右边的贡献，再做[mid+1,r]，这样可以保证前边的先完成，才能贡献后边。

• 我们考虑[l,mid]对于[mid+1,r]的贡献（这不就是CDQ分治吗？）。将[l,mid]的a提出来，与一个长度为r-l的b卷起来，得到答案再将对应位置贡献加到[mid+1,r]的对应位置中去。

• 每个位置会参与log次贡献，每贡献一次总共nlog(n)，均摊下来每个位置log的时间，所以总复杂度为O(nlog²n)

• 这题不失为一道分治NTT的模板题。

• 不会FFT/NTT的强烈推荐网上的一篇blog (实际上是我太菜了)

• https://www.cnblogs.com/zwfymqz/p/8244902.html

多项式求逆

• 将 Solution中的式子转化一下可以得到下列式子（注：将f[n]（上面是i）项右移，发现组合数和2的次幂这几个系数都为1，所以可以将sigma上的n-1变成n）
  
• 左式已知，右边仅有fk/(k-1)!那一项不知道，要求它。
• 设F=sigma(fk/(k-1)!),G=sigma(2^C(n-k,2)/(n-k)!)
• 考虑sigma可以看做卷积，卷积可以看成是多项式乘法。
• 所以上式可以写成H=F*G
• 已知多项式G，H，要求多项式F。
• 所以F=H*G^-1
• 将G多项式求逆一波后再用多项式乘法乘在一起即可。
• 多项式乘法也是一个NTT。
• 多项式求逆是一个倍增的NTT.
• 所以时间复杂度还是O(nlog²n)
• 有一个很好的多项式求逆总结的blog，挂在这里以供学习:
• https://www.cnblogs.com/yoyoball/p/8724115.html

分治NTT的代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#define maxn 130005
#define ll long long 
#define mo 1004535809
using namespace std;

ll n,i,j,k,f[maxn],g[maxn],_2[maxn],mul[maxn],invm[maxn];
ll a[maxn*2],b[maxn*2],bt[maxn*2];

ll ksm(ll x,ll y){
	ll s=1;
	for(;y;y/=2,x=x*x%mo) if (y&1)
		s=s*x%mo;
	return s;
}

void ntt(ll *a,int n,int sig){
	for(int i=0;i<n;i++) if (i<bt[i])
		swap(a[i],a[bt[i]]);
	for(int mid=1;mid<n;mid<<=1){
		ll gn=ksm(3,(mo-1)/(mid<<1));
		if (sig<0) gn=ksm(gn,mo-2);
		for(int j=0;j<n;j+=mid<<1){
			ll g=1;
			for(int k=0;k<mid;k++,g=g*gn%mo){
				ll x=a[j+k],y=g*a[j+k+mid]%mo;
				a[j+k]=(x+y)%mo;
				a[j+k+mid]=(x-y+mo)%mo;
			}
		}
	}
}

void merge(int l,int r){
	if (l==r) {f[l]=(f[l]+g[l])%mo;return;}
	int mid=(l+r)/2;
	merge(l,mid);
	int lim=1; while (lim<=r-l) lim<<=1;
	for(int i=1;i<lim;i++) bt[i]=(bt[i>>1]>>1)|((i&1)*(lim>>1));
	for(int i=0;i<lim;i++) a[i]=b[i]=0;
	for(int i=0;i<=mid-l;i++) a[i]=f[l+i]*invm[l+i-1]%mo;
	for(int i=0;i<=r-l;i++) b[i]=g[i]*invm[i]%mo;
	ntt(a,lim,1); ntt(b,lim,1);
	for(int i=0;i<lim;i++) a[i]=a[i]*b[i]%mo;
	ntt(a,lim,-1);
	ll ny=ksm(lim,mo-2);
	for(int i=mid+1;i<=r;i++) 
		f[i]=(f[i]-a[i-l]*ny%mo*mul[i-1]%mo+mo)%mo;
	merge(mid+1,r);
}

int main(){
	scanf("%lld",&n);
	_2[0]=1; for(i=1;i<=n;i++) _2[i]=_2[i-1]*2%mo;
	g[0]=1; for(i=1;i<=n;i++) g[i]=g[i-1]*_2[i-1]%mo;
	mul[0]=invm[0]=1;
	for(i=1;i<=n;i++) mul[i]=mul[i-1]*i%mo,invm[i]=invm[i-1]*ksm(i,mo-2)%mo;
	
	merge(1,n);
	printf("%lld",f[n]);
}

多项式求逆的代码

注意一个细节

• 当次数为lim的多项式A和B卷起来的时候，它们在做DFT和IDFT时的值域要开成lim * 2，因为它们的次数之和为lim * 2，尽管我们之后可能只需要次数lim以内的卷积结果，但是为了避免其他位置的印象要开大。
• 封装大法好。

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#define maxn 130005
#define maxm 520010
#define ll long long 
#define mo 1004535809
using namespace std;

ll n,i,j,k,lim,f[maxn],g[maxn],_2[maxn],mul[maxn],invm[maxn];
ll G[maxm],GG[maxm],H[maxm];
ll bt[maxm],A[maxm],B[maxm],C[maxm];

ll ksm(ll x,ll y){
	ll s=1;
	for(;y;y/=2,x=x*x%mo) if (y&1)
		s=s*x%mo;
	return s;
}

void ntt(ll *a,int n,int sig){
	for(int i=0;i<n;i++) if (i<bt[i])
		swap(a[i],a[bt[i]]);
	for(int mid=1;mid<n;mid<<=1){
		ll gn=ksm(3,(mo-1)/(mid<<1));
		if (sig<0) gn=ksm(gn,mo-2);
		for(int j=0;j<n;j+=mid<<1){
			ll g=1;
			for(int k=0;k<mid;k++,g=g*gn%mo){
				ll x=a[j+k],y=g*a[j+k+mid]%mo;
				a[j+k]=(x+y)%mo;
				a[j+k+mid]=(x-y+mo)%mo;
			}
		}
	}
}

ll a[maxm],b[maxm];

void mul_xl(ll *A,ll *B,ll lim){
	memset(a,0,sizeof(a));
	memset(b,0,sizeof(b));
	ll bat=lim<<1;
	for(int i=0;i<lim;i++) a[i]=A[i],b[i]=B[i];
	for(int i=0;i<bat;i++) bt[i]=(bt[i>>1]>>1)|((i&1)*(bat>>1));
	ntt(a,bat,1),ntt(b,bat,1);
	for(int i=0;i<bat;i++) a[i]=a[i]*b[i]%mo;
	ntt(a,bat,-1);
	ll ny=ksm(bat,mo-2);
	for(int i=0;i<bat;i++) A[i]=a[i]*ny%mo;
}

void inv(ll *A,ll *B,ll n){
	B[0]=ksm(A[0],mo-2);
	for(int lim=2;lim<=n*2;lim<<=1){
		memset(C,0,sizeof(C));
		for(int i=0;i<lim>>1;i++) C[i]=(B[i]<<1)%mo;
		mul_xl(B,B,lim>>1);
		mul_xl(B,A,lim);
		for(int i=0;i<lim;i++) B[i]=(C[i]-B[i]+mo)%mo;
	}
}

int main(){
	scanf("%lld",&n);
	_2[0]=1; for(i=1;i<=n;i++) _2[i]=_2[i-1]*2%mo;
	g[0]=1; for(i=1;i<=n;i++) g[i]=g[i-1]*_2[i-1]%mo;
	mul[0]=invm[0]=1;
	for(i=1;i<=n;i++) mul[i]=mul[i-1]*i%mo,invm[i]=invm[i-1]*ksm(i,mo-2)%mo;
	
	for(i=0;i<=n;i++) {
		G[i]=g[i]*invm[i]%mo;
		H[i]=g[i]*invm[i-1]%mo;
	}
	
	inv(G,GG,n);
	for(lim=1;lim<=n;lim<<=1);
	mul_xl(H,GG,lim);
	printf("%lld",H[n]*mul[n-1]%mo);
}