「bzoj省队十连测Day2 B」黑暗 - 图的计数+生成函数+分治FFT/多项式求逆

题目大意

    $n$个点的无向图，每条边都可能存在，一个图的权值是连通块个数的$m$次方，求所有可能的图的权值和。
    答案对$998244353$取模。
    多组数据。

---

数据范围

• 对于20%的数据$n\le7$
• 对于另外20%的数据$n\le500$
• 对于另外20%的数据$n\le2000$
• 对于100%的数据$T\le1000,n\le30000,m\le15$

---

题目分析

本题略难（KEKE秒了）。
我们对不同数据进行分析。

$20\%$算法

状压+并查集验证。

$40\%$算法

和城市规划类似。
设$g(i)$表示有$i$个点的有标号简单联通无向图个数，$h(i)$表示有$i$个点的有标号简单无向图个数。

$$h(n)=2^{\binom{n}2}$$

还是同样的强制$1$号点被包含，列出式子：

$$h(n)=\sum_{i=1}^n\binom{n-1}{i-1}g(i)h(n-i)$$

因此可以暴力计算$g(n)$。

设$f(i,j)$表示$i$个点$j$个连通块的个数。
枚举$1$号点所在连通块大小：

$$f(i,j)=\sum_{k=1}^if(i-k,j-1)g(k)\binom{i-1}{k-1}$$

时间复杂度为$O(n^3+nT)$。

$60\%$算法

可以直接设$f(i,j)$表示$i$个点的图，$j$个连通块，连通块数量的$j$次方之和。
仍然枚举$1$号点的连通块大小转移。
注意到：

$$(n+1)^m=\sum_{i=0}^m\binom min^i$$

因此得到转移式：

$$f(i,j)=\sum_{k=1}^{i}g(k)\binom{i-1}{k-1}\sum_{t=0}^j\binom jt f(i-k,t)$$

注意到$\sum_{t=0}^j\binom jt f(i-k,t)$与$i$没有直接关系，可以直接枚举$i-k,j$，在$O(nm^2)$的时间内预处理完毕。
总时间复杂度$O(n^2m+T)$。

$100\%$算法1

继续使用Stirling数展开自然数幂的套路。
将自然数幂转化为：

$$n^m=\sum_{k=0}^m{m \brace k}\binom nk k!$$

故可以转化为求$\binom{\text{连通块数}}m$。

假设$F(i,j)$表示$i$个点的图，权值是$\binom{\text{连通块数}}m$的权值和，同样枚举$1$所在连通块，注意到组合数的Pascal递推公式，可以得到：

$$F(i,j)=\sum_{k=1}^ig(k)\binom{i-1}{k-1}(F(i-k,j)+F(i-k,j-1))$$

则：

$$f(n,m)=\sum_{k=0}^m{m \brace k}F(n,k)k!$$

Stirling数很好预处理，现在考虑如何快速求出$F$。
将组合数展开，得到：

$$\frac{F(i,j)}{(i-1)!}=\sum_{k=1}^i\frac{g(k)}{(k-1)!}\frac{F(i-k,j)+F(i-k,j-1)}{(i-k)!}$$

如何就可以使用分治FFT计算了。
时间复杂度$O(mn\log^2 n+mT)$。

$100\%$算法2

官方题解看不懂，因为图的计数包含标号重分配的过程，因此使用指数型生成函数推导。
令指数型生成函数$F(i)$表示$i$个点连通图个数，$G(i)$表示$i$个点任意图个数。
那么要求的答案即为：

$$\sum_{i=0}^{+\infty}\frac{i^mF^i}{i!}$$

而根据麦克劳林级数有：

$$e^x=\sum_{i=0}^{+\infty}\frac{x^i}{i!}$$

因此可得：

$$e^F=\sum_{i=0}^{+\infty}\frac{F^i}{i!}=G$$

前方公式恐惧症慎入

$$\begin{aligned} ans&=\sum_{i=0}^{+\infty}\frac{i^mF^i}{i!} \\ &=\sum_{i=0}^{+\infty}\frac{\sum\limits_{k=0}^m{m \brace k}\binom ikk!F^i}{i!} \\ &=\sum_{i=0}^{+\infty}\frac{\sum\limits_{k=0}^m{m \brace k}\frac{i!}{(i-k)!}F^i}{i!} \\ &=\sum_{i=0}^{+\infty}(\sum\limits_{k=0}^m{m \brace k}\frac{F^i}{(i-k)!}) \\ &=\sum_{i=0}^{+\infty}(\sum\limits_{k=0}^m{m \brace k}F^k\frac{F^{i-k}}{(i-k)!}) \\ &=\sum_{k=0}^m{m \brace k}F^k\sum_{i=0}^{+\infty}\frac{F^{i-k}}{(i-k)!} \\ &=\sum_{k=0}^m{m \brace k}F^kG \end{aligned}$$

$m=15$，$k$次幂可以暴力预处理，故瓶颈在于求出$F$。

因为$e^F=G$，故$F=\ln G$。
$G$是很好求出的，故现在只需要对$G$做一个多项式$lnp$即可。
利用求导再积分的方法可以算出多项式$lnp$，因此问题可以很简单的解决，时间复杂度$O(nm\log n+T)$。
因为常数原因，实际比分治FFT跑得慢。

---

代码

#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;

inline const int Get_Int() {
	int num=0,bj=1;
	char x=getchar();
	while(x<'0'||x>'9') {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(x>='0'&&x<='9') {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}

const int maxn=65536+5,maxk=17;
const LL mod=998244353,g=3;

LL Quick_Pow(LL a,LL b) {
	LL sum=1;
	for(; b; b>>=1,a=a*a%mod)if(b&1)sum=sum*a%mod;
	return sum;
}

LL inv(LL x) {
	return Quick_Pow(x,mod-2);
}

struct NumberTheoreticTransform {
	int n,rev[maxn];
	LL omega[maxn],iomega[maxn];
	
	void init(int n) {
		this->n=n;
		int x=Quick_Pow(g,(mod-1)/n);
		omega[0]=iomega[0]=1;
		for(int i=1; i<n; i++) {
			omega[i]=omega[i-1]*x%mod;
			iomega[i]=inv(omega[i]);
		}
		int k=log2(n);
		for(int i=0; i<n; i++) {
			int t=0;
			for(int j=0; j<k; j++)if(i&(1<<j))t|=(1<<(k-j-1));
			rev[i]=t;
		}
	}
	
	void transform(LL* a,LL* omega) {
		for(int i=0; i<n; i++)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int len=2; len<=n; len*=2) {
			int mid=len>>1;
			for(LL* p=a; p!=a+n; p+=len)
				for(int i=0; i<mid; i++) {
					LL t=omega[n/len*i]*p[mid+i]%mod;
					p[mid+i]=(p[i]-t+mod)%mod;
					p[i]=(p[i]+t)%mod;
				}
		}
	}

	void dft(LL* a) {
		transform(a,omega);
	}

	void idft(LL* a) {
		transform(a,iomega);
		LL x=inv(n);
		for(int i=0; i<n; i++)a[i]=a[i]*x%mod;
	}
} ntt;

void polynomial_inverse(const LL* a,const int n,LL* b) {
	if(n==1) {
		b[0]=inv(a[0]);
		return;
	}
	polynomial_inverse(a,n>>1,b);
	int p=n<<1;
	static LL x[maxn];
	copy(a,a+n,x),fill(x+n,x+p,0);
	ntt.init(p),ntt.dft(x),ntt.dft(b);
	for(int i=0; i<p; i++)b[i]=b[i]*((2-x[i]*b[i]%mod+mod)%mod)%mod;
	ntt.idft(b),fill(b+n,b+p,0);
}

void Multiply(const LL* a1,const int n1,const LL* a2,const int n2,LL* ans) {
	int n=1;
	while(n<n1+n2)n<<=1;
	ntt.init(n);
	static LL c1[maxn],c2[maxn];
	copy(a1,a1+n1,c1),fill(c1+n1,c1+n,0);
	copy(a2,a2+n2,c2),fill(c2+n2,c2+n,0);
	ntt.dft(c1);
	ntt.dft(c2);
	for(int i=0; i<n; i++)c1[i]=c1[i]*c2[i]%mod;
	ntt.idft(c1);
	for(int i=0; i<n1+n2-1; i++)ans[i]=c1[i];
}

LL tmp[maxn];

void polynomial_lnp(LL *a,LL *ans,int n) {
	int p=1;
	while(p<n)p<<=1;
	polynomial_inverse(a,p,tmp);
	fill(tmp+n,tmp+p,0);
	for(int i=0; i<n; i++)ans[i]=a[i+1]*(i+1)%mod; //求导
	ans[n]=0;
	Multiply(ans,n,tmp,n,ans);
	fill(ans+n,ans+2*n,0);
	for(int i=n; i>=1; i--)ans[i]=ans[i-1]*inv(i)%mod; //积分
	ans[0]=0;
}

int n;
LL S[maxk][maxk],fac[maxn],invf[maxn],G[maxn],F[maxn],Mul[maxk][maxn];

int main() {
	n=30000;
	S[0][0]=1;
	for(int i=1; i<maxk; i++) {
		for(int j=1; j<maxk; j++)S[i][j]=(S[i-1][j-1]+S[i-1][j]*j)%mod;
	}
	fac[0]=invf[0]=1;
	for(int i=1; i<=n; i++) {
		fac[i]=fac[i-1]*i%mod;
		invf[i]=inv(fac[i]);
	}
	for(int i=0; i<=n; i++)G[i]=Quick_Pow(2,i*(i-1)>>1)*invf[i]%mod;
	polynomial_lnp(G,F,n+1);
	Mul[0][0]=1;
	for(int i=1; i<maxk; i++) {
		Multiply(F,n+1,Mul[i-1],n+1,Mul[i]);
		fill(Mul[i]+n+1,Mul[i]+2*n,0);
	}
	for(int i=maxk-1; i>=1; i--) {
		for(int j=1; j<i; j++)
			for(int k=0; k<=n; k++)Mul[i][k]=(Mul[i][k]+Mul[j][k]*S[i][j])%mod;
		Multiply(Mul[i],n+1,G,n+1,Mul[i]);
		fill(Mul[i]+n+1,Mul[i]+2*n,0);
	}
	int t=Get_Int();
	while(t--) {
		int n=Get_Int(),m=Get_Int();
		printf("%lld\n",Mul[m][n]*fac[n]%mod);
	}
	return 0;
}