【UOJ #50】【UR #3】—链式反应（生成函数+分治NTT/多项式Exp+常微分方程）

传送门

---

首先这个题面就很胃疼
而且感觉讲的不是很清楚

实际上是要求满足如下条件的树的个数：

对于每个非叶节点，有$c+2$个儿子，其中有$cin A$个叶子节点和2个非叶节点
且点的编号满足父亲小于儿子

按照套路设$f[i]$为$i$个点的答案可以列出$dp$式

$f[i]=frac 1 2sum_{j}sum_{k}[i-j-k-1in A]{i-1choose j}{i-j-1choose k}f_jf_k$
中间$frac 1 2$是由于两个非叶儿子会把每种情况算2次
边界$f[1]=1$

设$A(x)=sum_{iin A}frac{x^i}{i!},f(x)=sum_{i}f[i]frac{x^i}{i!}$

那么有$sum_{i}f[i+1]frac{x^i}{i+1}=f'=frac 1 2f^2A+1$

这个可以分治$NTT$做到$O(nlog^2n)$而且跑的比$nlogn$快

不过注意如果$l ot=1$的话实际上情况只统计到了一次
而前面有一个$frac 1 2$,所以算的时候需要乘2

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define pb push_back
#define re register
#define ll long long
#define pii pair<int,int>
#define fi first
#define bg begin
#define se second
#define poly vector<int>
cs int RLEN=1<<20|1;
inline char gc(){
	static char ibuf[RLEN],*ib,*ob;
	(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
	return (ib==ob)?EOF:*ib++;
}
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0;bool f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=998244353,G=3;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return r>=mod?r%mod:r;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=r>=mod?r%mod:r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
cs int C=19;
int *w[C+1],rev[(1<<C)|5];
inline void init_w(){
	for(int i=1;i<=C;i++)w[i]=new int[1<<(i-1)];
	int wn=ksm(G,(mod-1)/(1<<C));
	w[C][0]=1;
	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--)
	for(int j=0;j<(1<<(i-1));j++)w[i][j]=w[i+1][j<<1];
}
inline void init_rev(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1,a0,a1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=mid<<1)
	for(int j=0;j<mid;j++)
	a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1){
		reverse(f.bg()+1,f.bg()+lim);
		for(int i=0,iv=Inv(lim);i<lim;i++)Mul(f[i],iv);
	}
}
inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=32){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
		Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;
	init_rev(lim);
	a.resize(lim),b.resize(lim);
	ntt(a,lim,1),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}
cs int N=200005;
int inv[N<<1],fac[N<<1],ifac[N<<1];
inline void init_inv(){
	cs int len=(N-5)<<1;
	fac[0]=ifac[0]=inv[0]=inv[1]=1;
	for(int i=2;i<=len;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
	for(int i=1;i<=len;i++)fac[i]=mul(fac[i-1],i),ifac[i]=mul(ifac[i-1],inv[i]);
}
int a[N],g[N],f[N],n;
char s[N];
void cdq(int l,int r){
	if(l==r)return;
	int mid=(l+r)>>1,l1=r-l+1;
	cdq(l,mid);
	poly h,p;
	for(int i=l;i<=mid;i++)h.pb(f[i]);
	for(int i=1;i<=l1;i++)p.pb(f[i]);
	h=h*p;int coef=1+(l>1);
	for(int i=mid+1;i<=r;i++)Add(g[i],mul(coef,h[i-l-1]));
	h.clear(),p.clear();
	for(int i=l;i<=mid;i++)h.pb(g[i]);
	for(int i=1;i<=l1;i++)p.pb(a[i]);
	h=h*p; 
	for(int i=mid+1;i<=r;i++)Add(f[i],mul(inv[i*2],h[i-l-1]));
	cdq(mid+1,r);
}
int main(){
	#ifdef Stargazer
	freopen("lx.cpp","r",stdin);
	#endif
	init_w(),init_inv();
	n=read();
	scanf("%s",s+1);
	for(int i=1;i<=n;i++)a[i]=mul((s[i]=='1'),ifac[i-1]);
	f[1]=1,cdq(1,n);
	for(int i=1;i<=n;i++)cout<<mul(f[i],fac[i])<<'
';
}

---

接下来是$O(nlogn)$的做法
而且是解常微分方程的一个通用套路

先设$h(x)=frac 1 2Ax^2+1$
那么实际要解的就是$f'(x)=h(f(x))$
按照牛顿迭代的通用方法
假设已经求出$f_0'(x)equiv h(f_0(x))\%x^{lceil frac x 2 ceil}$
求$f'(x)equiv h(f(x))\% x^n$
那么$f'(x)=h(f_0)+h'(f_0)(f-f_0)$
由于这个东西无法方便的牛顿迭代
考虑构造$v(x)=exp(-int h'(f_0(x))mathrm{dx})$
那么有$v'=-h'(f_0)v$
等式两边乘上$v$
有$f'v=v(h(f_0)+h'(f_0)(f-f_0))$
$f'v-fh'(f_0)v=v(h(f_0)-h'(f_0)f_0)$
$f'v+fv'=v(h(f_0)-h'(f_0)f_0)$
$(fv)'=v(h(f_0)-h'(f_0)f_0)$
$f=frac{1}{v}int(v(h(f_0)-h'(f_0)f_0))mathrm{dx}$

因为这个和具体$h$没有关系所以可以通用

对于原来的$h$

那么就有$f=frac 1 vint(v(1-frac 1 2Af_0^2))mathrm{dx}$

复杂度$O(nlogn)$跑的远没有$nlog^2n$快

这东西知道一下就好了
除了理论复杂度以外其他都比不过分治$ntt$

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define pb push_back
#define re register
#define ll long long
#define pii pair<int,int>
#define fi first
#define bg begin
#define se second
#define poly vector<int>
cs int RLEN=1<<20|1;
inline char gc(){
	static char ibuf[RLEN],*ib,*ob;
	(ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
	return (ib==ob)?EOF:*ib++;
}
#define gc getchar
inline int read(){
	char ch=gc();
	int res=0;bool f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
	return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=998244353,G=3;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return r>=mod?r%mod:r;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=r>=mod?r%mod:r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
cs int C=20;
int *w[C+1],rev[(1<<C)|5],inv[(1<<C)+5];
inline void init_w(){
	for(int i=1;i<=C;i++)w[i]=new int[1<<(i-1)];
	int wn=ksm(G,(mod-1)/(1<<C));
	inv[0]=inv[1]=w[C][0]=1;
	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--)
	for(int j=0;j<(1<<(i-1));j++)w[i][j]=w[i+1][j<<1];
	for(int i=2,lim=(1<<(C-1));i<lim;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
}
inline void init_rev(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1,a0,a1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=mid<<1)
	for(int j=0;j<mid;j++)
	a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1){
		reverse(f.bg()+1,f.bg()+lim);
		for(int i=0,iv=Inv(lim);i<lim;i++)Mul(f[i],iv);
	}
}
inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=32){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
		Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;
	init_rev(lim);
	a.resize(lim),b.resize(lim);
	ntt(a,lim,1),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}
inline poly Inv(poly a,int deg){
	poly b(1,Inv(a[0])),c;
	for(int lim=4;lim<(deg<<2);lim<<=1){
		c.resize(lim>>1);
		for(int i=0;i<(lim>>1);i++)c[i]=(i<a.size())?a[i]:0;
		init_rev(lim);
		b.resize(lim),c.resize(lim);
		ntt(c,lim,1),ntt(b,lim,1);
		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(c[i],b[i])));
		ntt(b,lim,-1),b.resize(lim>>1);
	}b.resize(deg);return b;
}
inline poly deriv(poly a){
	for(int i=0;i+1<a.size();i++)a[i]=mul(a[i+1],i+1);
	a.pop_back();return a;
}
inline poly integ(poly a){
	a.pb(0);
	for(int i=(int)a.size()-1;~i;i--)a[i]=mul(a[i-1],inv[i]);
	a[0]=0;return a;
}
inline poly Ln(poly a,int deg){
	a=integ(Inv(a,deg)*deriv(a)),a.resize(deg);return a;
}
inline poly exp(poly a,int deg){
	poly b(1,1),c;
	for(int lim=2;lim<(deg<<1);lim<<=1){
		c=Ln(b,lim);
		for(int i=0;i<lim;i++)c[i]=dec((i<a.size())?a[i]:0,c[i]);
		Add(c[0],1);
		b=b*c,b.resize(lim);
	}b.resize(deg);
	return b;
}
cs int N=200010;
int n,a[N],fac[N],ifac[N];
inline poly calc(int lim){
	if(lim==1)return poly(1,0);
	poly f0=calc((lim+1)>>1);
	poly p,v;
	for(int i=0;i<lim;i++)p.pb(a[i]);
	p=p*f0,v=integ(p),v.resize(lim);
	for(int i=0;i<lim;i++)v[i]=dec(0,v[i]);
	v=exp(v,lim);
	poly h=p*f0;h.resize(lim);
	for(int i=0;i<lim;i++)h[i]=dec(0,mul(inv[2],h[i]));
	Add(h[0],1),h=integ(v*h);
	h=h*Inv(v,lim),h.resize(lim);
	return h;
}
inline void init_inv(){
	cs int len=N-5;
	fac[0]=ifac[0]=1;
	for(int i=1;i<=len;i++)fac[i]=mul(fac[i-1],i);
	ifac[len]=Inv(fac[len]);
	for(int i=len-1;i;i--)ifac[i]=mul(ifac[i+1],i+1);
}
char s[N];
int main(){
	#ifdef Stargazer
	freopen("lx.cpp","r",stdin);
	#endif
	init_w(),init_inv();
	n=read();
	scanf("%s",s);
	for(int i=0;i<n;i++)a[i]=mul(ifac[i],s[i]=='1');
	poly ans=calc(n+1);
	for(int i=1;i<=n;i++)cout<<mul(fac[i],ans[i])<<'
';
	return 0;
}