BZOJ 3456 城市规划 (组合计数、DP、FFT)

题目链接: https://www.lydsy.com/JudgeOnline/problem.php?id=3456

著名的多项式练习题，做法也很多，终于切掉了纪念

首先求一波递推式: 令$F(n)$为$n$个点的有标号无向连通图的个数，则考虑补集转化为有标号无向不连通图的个数，然后枚举$1$号点所在联通块的大小: $$F(n)=2^{n\choose 2}-\sum^{n-1}_{i=1} {n-1\choose i-1} F(i)2^{n-i\choose 2}$$

这样可以做到$O(n^2)$, 后面就该大佬们各显神通了，我在这里整理一下四种做法:

做法一

直接使用分治NTT优化，时间复杂度$O(n\log^2n)$。但我不会分治NTT，所以不具体说了。

做法二

\[F(n)=2^{n\choose 2}-\sum^{n-1}_{i=1}\frac{(n-1)!}{(i-1)!(n-i)!}F(i)2^{n-i\choose 2}\\
\frac{F(n)}{(n-1)!}=\frac{2^{n\choose 2}}{(n-1)!}-\sum^{n-1}_{i=1} \frac{F(i)}{(i-1)!}\frac{2^{n-i\choose 2}}{(n-i)!}$$移项可得$$\frac{2^{n\choose 2}}{(n-1)!}=\sum^{n}_{i=1} \frac{F(i)}{(i-1)!}\frac{2^{n-i\choose 2}}{(n-i)!}\]

令$A(x)=\sum_{n>0}\frac{F(n)}{(n-1)!}, G(x)=\sum_{n\ge 0}\frac{2^{n\choose 2}}{n!}, H(x)=\sum_{n>0}\frac{2^{n\choose 2}}{(n-1)!}$, 则有$$F(x)G(x)=H(x)\ F(x)=H(x)G(x)^{-1}$$

多项式求逆即可。

时间复杂度$O(n\log n)$.

这应该是代码复杂度和运行效率上最好的一种做法，但是做法三和做法四也有一定的启发意义。

做法三

设$G(n)=2^{n\choose 2}$表示$n$个点有标号无向图的个数。设$F(x),G(x)$分别为$F(n),G(n)$的指数生成函数(EGF).

由于一个有标号无向图由若干个彼此之间无顺序的联通块组成，因此其指数生成函数$G(x)=\sum_{n\ge 1}\frac{F(x)^n}{n!}$.

即$G(x)=e^{F(x)}$, $F(x)=\ln G(x)$. 多项式$\ln$即可。

时间复杂度$O(n\log n)$.

做法四

(这个做法是我自己想的，有错敬请指出)(这种做法其实是用另一种方式推导做法三)

感谢_rqy大爷博客里的生成函数简介。

仿照求Bell数的EGF方法，进行以下推导: $$\frac{F(n)}{n!}=\frac{G(n)}{n!}-\sum^{n-1}{i=1} \frac{F(i)}{n(i-1)!}\frac{G(n-i)}{(n-i)!}\ \frac{G(n)}{n!}=\frac{F(n)}{n!}+\frac{1}{n}\sum^{n-1}{i=1}\frac{iF(i)}{i!}\frac{G(n-i)}{(n-i)!}$$

这里我们发现$\frac{iF(i)}{i!}$就是$[x^{i-1}]F'(x)$, 于是上式可以改写为$$[x^n]G(x)=[xn]F(x)+\frac{1}{n}\sum^{{n-1}_{i=1}[x}{i-1}]F'(x)\times [x^{n-i}]G(x)\ =\frac{1}{n}(n[x^n]F(x)+\sum{n-1}_{i=1}[x^{i-1}]F'(x)\times [x^{n-i}]G(x))\ =\frac{1}{n}\sum^{n}_{i=1}[x{i-1}]F'(x)\times [x^{n-i}]G(x)\ G(x)=\int^x_0 F'(x)G(x)\text{d}x\ \frac{G'(x)}{G(x)}=F'(x)\ \ln G(x)=F(x)$$.

一定要注意求和边界! 我推式子的时候没注意求和上界是$n$还是$(n-1)$的问题结果一直推出来$G(x)=F(x)+\int^x_0 F'(x)G(x)\text{d}x$查了一小时……

代码

做法二

#include<cstdio>

#include<cstdlib>

#include<cstring>

#include<cassert>

#include<iostream>

#define llong long long

using namespace std;

inline int read()

{

	int x=0; bool f=1; char c=getchar();

	for(;!isdigit(c);c=getchar()) if(c=='-') f=0;

	for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0');

	if(f) return x;

	return -x;

}

const int N = 1<<19;

const int LGN = 19;

const int P = 1004535809;

const int G = 3;

llong fact[N+3],finv[N+3];

llong quickpow(llong x,llong y)

{

	llong cur = x,ret = 1ll;

	for(int i=0; y; i++)

	{

		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}

		cur = cur*cur%P;

	}

	return ret;

}

llong mulinv(llong x) {quickpow(x,P-2);}

namespace Polynomial

{

	llong tmp1[N+3],tmp2[N+3],tmp3[N+3],tmp4[N+3],tmp5[N+3],tmp6[N+3];

	int fftid[N+3];

	int getdgr(int n)

	{

		int ret = 1; while(ret<=n) ret<<=1;

		return ret;

	}

	void init_fftid(int dgr)

	{

		int len = 0; for(int i=1; i<=LGN; i++) {if((1<<i)==dgr) {len = i; break;}}

		for(int i=1; i<dgr; i++) fftid[i] = (fftid[i>>1]>>1)|((i&1)<<(len-1));

	}

	void ntt(int dgr,int coe,llong poly[],llong ret[])

	{

		init_fftid(dgr);

		if(poly==ret) {for(int i=0; i<dgr; i++) {if(i<fftid[i]) swap(ret[i],ret[fftid[i]]);}}

		else {for(int i=0; i<dgr; i++) ret[i] = poly[fftid[i]];}

		for(int i=1; i<=(dgr>>1); i<<=1)

		{

			llong tmp = quickpow(G,(P-1)/(i<<1));

			if(coe==-1) tmp = mulinv(tmp);

			for(int j=0; j<dgr; j+=(i<<1))

			{

				llong expn = 1ll;

				for(int k=0; k<i; k++)

				{

					llong x = ret[j+k],y = ret[i+j+k]*expn%P;

					ret[j+k] = (x+y)%P;

					ret[j+i+k] = (x-y+P)%P;

					expn = expn*tmp%P;

				}

			}

		}

		if(coe==-1)

		{

			llong tmp = mulinv(dgr);

			for(int i=0; i<dgr; i++) ret[i] = ret[i]*tmp%P;

		}

	}

	void polymul(int dgr,llong poly1[],llong poly2[],llong ret[])

	{

		ntt(dgr<<1,1,poly1,tmp1); ntt(dgr<<1,1,poly2,tmp2);

		for(int i=0; i<(dgr<<1); i++) tmp2[i] = tmp1[i]*tmp2[i]%P;

		ntt(dgr<<1,-1,tmp2,ret);

	}

	void polyinv(int dgr,llong poly[],llong ret[])

	{

		for(int i=0; i<dgr; i++) ret[i] = tmp1[i] = 0ll;

		ret[0] = mulinv(poly[0]); tmp1[0] = poly[0];

		for(int i=1; i<=(dgr>>1); i<<=1)

		{

			for(int j=i; j<(i<<1); j++) tmp1[j] = poly[j];

			ntt((i<<2),1,ret,tmp2); ntt((i<<2),1,tmp1,tmp3);

			for(int j=0; j<(i<<2); j++) tmp2[j] = tmp2[j]*tmp2[j]%P*tmp3[j]%P;

			ntt((i<<2),-1,tmp2,tmp3);

			for(int j=0; j<(i<<1); j++) ret[j] = (2ll*ret[j]-tmp3[j]+P)%P;

		}

	}

}

llong f[N+3],g[N+3],h[N+3],ginv[N+3];

int n;

int main()

{

	fact[0] = 1ll; for(int i=1; i<=N; i++) fact[i] = fact[i-1]*i%P;

	finv[N] = quickpow(fact[N],P-2); for(int i=N-1; i>=0; i--) finv[i] = finv[i+1]*(i+1)%P;

	scanf("%d",&n); int dgr = Polynomial::getdgr(n);

	for(int i=0; i<=n; i++) {g[i] = quickpow(2ll,i*(i-1ll)/2ll)*finv[i]%P;}

	for(int i=1; i<=n; i++) {h[i] = quickpow(2ll,i*(i-1ll)/2ll)*finv[i-1]%P;}

//	printf("g: "); for(int i=0; i<dgr; i++) printf("%lld ",g[i]); puts("");

//	printf("h: "); for(int i=0; i<dgr; i++) printf("%lld ",h[i]); puts("");

	Polynomial::polyinv(dgr,g,ginv);

//	printf("ginv: "); for(int i=0; i<dgr; i++) printf("%lld ",ginv[i]); puts("");

	Polynomial::polymul(dgr,ginv,h,f);

	printf("%lld\n",f[n]*fact[n-1]%P);

	return 0;

}

生成函数这东西真的是有趣！！！

巴特西

BZOJ 3456 城市规划 (组合计数、DP、FFT)

做法一

做法二

做法三

做法四

代码

最新文章

热门文章