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

最新推荐文章于 2025-01-22 09:24:34 发布

YiPeng_Deng

最新推荐文章于 2025-01-22 09:24:34 发布

阅读量427

点赞数 1

分类专栏：题解多项式 FFT和NTT 组合数学容斥计数学习小计文章标签：分治NTT 多项式求逆

本文链接：https://blog.csdn.net/qq_43649416/article/details/94867106

版权

题解同时被 3 个专栏收录

133 篇文章

订阅专栏

学习小计

34 篇文章

订阅专栏

计数

32 篇文章

订阅专栏

本文探讨了如何利用图论中的连通图概念和算法解决城市规划问题，具体包括使用分治NTT和多项式求逆等高级算法来计算特定规模的无向连通图数目。文章详细解释了算法原理，并提供了具体的代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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);
}