BZOJ 4555: [Tjoi2016&Heoi2016]求和 [分治FFT 组合计数 | 多项式求逆]

4555: [Tjoi2016&Heoi2016]求和

题意：求$$
sum_{i=0}^n sum_{j=0}^i S(i,j)cdot 2^jcdot j!
S是第二类斯特林数

[ *** 首先你要把这个组合计数肝出来，~~于是我去翻了一波《组合数学》~~ *用斯特林数容斥原理推导那个式子可以直接出卷积形式，见下一篇，本篇是分治fft做法* </br> ###组合计数 **斯特林数** $S(n,i)$表示将n个不同元素划分成i个相同集合非空的方案数 **Bell数** $B(n)=sumlimits_{i=0}^n S(n,i)$也就是分成任意个相同的集合 有一个递推式 ]

B(n) = sum_{i=0}^{n-1} inom{n-1}{n-1-i}B(i) = sum_{i=0}^{n-1} inom{n-1}{i}B(i)

[$推导$： 考虑和第一个元素在同一个集合的元素有几个，剩下的是子问题 或者这么想，因为每个集合相同，我们规定第一个元素在一个集合后这个集合就不相同了，剩下的元素分入这个不相同的集合或者成为子问题 </br> $S'(n,i)=S(n,i)*i!$就是**集合不相同**啦， ]

B'(n) = sum_{i=1}^n inom{n}{i}B'(n-i)

[$推导$： 因为每个集合都不相同，我们没有必要规定第一个元素在集合里了，直接枚举第一个集合元素个数就可以了 </br> 然后那个$2^j$可以认为每个集合有两种颜色，枚举当前集合顺便枚举颜色，乘上一个2就行了 令$f(n)=2 sumlimits_{i=1}^n inom{n}{i} f(n-i)$，答案就是$sum_{i=0}^n f(i)$ 快速求出f就可以做了，整理一下发现 ]

f(n) = 2 n! sum_{i=1}^n frac{1}{i!} frac{f(n-i)}{(n-i)!}

[是一个卷积的形式，但是两边都有f，所以要用**CDQ分治套FFT** </br> ###分治fft 和CDQ分治一样，每次用$[l,mid]$更新$[mid+1,r]$ 可以发现，$f(i)[mid+1,r]$都是$frac{f(i)}{i!}[l,mid]$卷上$frac{1}{i!}[1,r-l]$的结果 我们把两个函数$[l,mid]$和$[1,r-l]$的部分拿出来做卷积，更新$[mid+1,r]$ 可以认为先把向量移动了l位 复杂度$O(nlog^2n)$，13560ms，~~貌似别人的分治fft更慢20000ms多~~ </br> 注意初始值$f[0]=1$因为$S(0,0)=1$，以及cdq(0,n) </br> ###多项式求逆 当然也可以，具体看[51nod这篇讨论](http://www.51nod.com/question/index.html#!questionId=1713)吧 把$f_0 x^0$单独拿出来的思想比较有意思,$A(x)$从0开始，$B(x)$从1开始（$B_0 = 0$） $$A(x) = F_0 x^0 + A(x) B(x)]

update 2017.5.3 : 发现tls貌似补了一句(2x e^x)，求这个东西来做也可以，但不如直接处理处B来方便

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=(1<<18)+5, INF=1e9;
const ll P=998244353, g=3;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}

ll Pow(ll a, ll b) {
	ll ans=1;
	for(; b; b>>=1, a=a*a%P)
		if(b&1) ans=ans*a%P;
	return ans;
}

int n, rev[N];
ll inv[N], fac[N], facInv[N];
ll f[N], a[N], b[N];
void dft(ll *a, int n, int flag) {
	for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
	for(int l=2; l<=n; l<<=1) {
		int m=l>>1;
		ll wn = Pow(g, flag==1 ? (P-1)/l : P-1-(P-1)/l);
		for(ll *p=a; p!=a+n; p+=l) {
			ll w=1;
			for(int k=0; k<m; k++) {
				ll t = w*p[k+m];
				p[k+m]=(p[k]-t+P)%P;
				p[k]=(p[k]+t)%P;
				w=w*wn%P;
			}
		}
	}
	if(flag==-1) {
		ll inv=Pow(n, P-2);
		for(int i=0; i<n; i++) a[i]=a[i]*inv%P;
	}
}
void cdq(int l, int r) {
	if(l==r) return;
	int mid=(l+r)>>1;
	cdq(l, mid);
	int lim=r-l+1, n=1, k=0;
	while(n<lim) n<<=1, k++;
	for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));

	for(int i=0; i<n; i++) a[i]=b[i]=0;
	for(int i=l; i<=mid; i++) a[i-l] = f[i];
	for(int i=0; i<=r-l; i++) b[i] = facInv[i];
	dft(a, n, 1); dft(b, n, 1);
	for(int i=0; i<n; i++) a[i]=a[i]*b[i]%P;
	dft(a, n, -1);

	for(int i=mid+1; i<=r; i++) f[i]=(f[i] + 2*a[i-l])%P;
	cdq(mid+1, r);
}
int main() {
	freopen("in","r",stdin);
	n=read();
	inv[1]=1; fac[0]=facInv[0]=1;
	for(int i=1; i<=n; i++) {
		if(i!=1) inv[i] = (P-P/i)*inv[P%i]%P;
		fac[i] = fac[i-1]*i%P;
		facInv[i] = facInv[i-1]*inv[i]%P;
	}
	f[0]=1;
	cdq(0, n);
	ll ans=0;
	for(int i=0; i<=n; i++) ans = (ans + f[i]*fac[i]%P)%P;
	if(ans<0) ans+=P;
	printf("%lld
", ans);
}