【集训队作业2018】普通的计数题（牛顿迭代）（微分方程）（多项式全家桶）

传送门

• 首先完成模型转换：一个 1 操作的点向消去的连边，那么有可能连的全部是 0 边，或是有一些 1 边和 0 边，发现这个就是一棵树，结点满足大的在上，如果一个点的儿子全是叶子，那么儿子个数 $\in B$，否则叶子儿子个数 $\in A$
• 我们直接暴力上 $EGF$，树的 $EGF$ 记为 $F(x)$，集合的 $EGF$ 记为 $A,B$
  为了方便我们强行往 $B$ 中填一个 0，那么显然有 $f_0=0,f_1=1$
  $$\begin{gathered} \frac{f_n}{(n-1)!}=\frac{[n-1\in B]}{(n-1)!}\sum _{i\in A}\frac{1}{i!}\sum _{k\ge 1}\frac{1}{k!}\sum _{{\sum }_{i=1}^k{a}_j=n-i-1,a_j\ge 2}\prod \frac{f_{a_j}}{a_j!} \\ \Rightarrow F^{\prime }=B+(A\exp (F-x)-1),F^{\prime }=A{e}^{-x}\exp F+(B-A) \end{gathered}$$
  即解一个微分方程 $F^{\prime }\equiv h(F)(\text{mod}{x}^n)$，假设已经求得 $F_0$ 使得
  $F_0^{\prime }\equiv h(F_0)(\text{mod}{x}^{\lceil \frac{x}{2}\rceil })$，泰勒展开，$F^{\prime }=h(F_0)+h^{\prime }(F_0)(F-F_0)(\text{mod}{x}^n)$
  有 $F^{\prime }-h^{\prime }(F_0)F=h(F_0)-h^{\prime }(F_0)F_0$，注意到我们可以构造一个函数 $G$ 使得 $(FG{)}^{\prime }=F^{\prime }G+G^{\prime }F=F^{\prime }G+(-h^{\prime }(F_0)G)F$，那么 $G=\exp (-\int h^{\prime }(F_0)\text{d}x)$
  于是我们将方程两边乘上 $G$ 再积分得
  $$F=\frac{1}{G}\int G(h(F_0)-h^{\prime }(F_0)F_0)\text{d}x$$
  这样就可以迭代了，复杂度 $O(nlog(n))$

#include<bits/stdc++.h>
#define cs const
#define poly vector<int>
#define pb push_back
using namespace std;
int read(){
	int cnt=0, f=1; char ch=0;
	while(!isdigit(ch)){ ch=getchar(); if(ch=='-') f=-1; }
	while(isdigit(ch)) cnt=cnt*10+(ch-'0'), ch=getchar();
	return cnt*f;
}
cs int Mod = 998244353;
int add(int a, int b){ return a + b >= Mod ? a + b - Mod : a + b; }
int dec(int a, int b){ return a - b < 0 ? a - b + Mod : a - b; }
int mul(int a, int b){ return 1ll * a * b % Mod; }
int ksm(int a, int b){ int as=1; for(;b;b>>=1,a=mul(a,a)) if(b&1) as=mul(as,a); return as; }
void Add(int &a, int b){ a = add(a, b); }
void Mul(int &a, int b){ a = mul(a, b); }
void Dec(int &a, int b){ a = dec(a, b); }
cs int N = 2e5 + 50;
cs int M = 1 << 18 | 5;
cs int K = 18;
int n, ifac[N], fac[N], iv[M];
void pre_work(int n){
	iv[0]=iv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1; 
	for(int i=2; i<=(1<<K); i++) iv[i]=mul(Mod-Mod/i,iv[Mod%i]);
	for(int i=2; i<=n; i++) fac[i]=mul(fac[i-1],i), ifac[i]=mul(ifac[i-1],iv[i]);
}
poly w[K+1];
int up, bit; poly rev;
void output(poly a){
	for(int i=0; i<a.size(); i++) cout<<a[i]<<" "; puts("");
}
void init(int deg){
	up=1, bit=0; while(up<deg) up<<=1,++bit; rev.resize(up); 
	for(int i=0; i<up; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void NTT_init(){
	for(int i=1; i<=K; i++) w[i].resize(1<<(i-1));
	int wn=ksm(3,(Mod-1)/(1<<K)); w[K][0]=1;
	for(int i=1; i<(1<<(K-1)); i++) w[K][i]=mul(w[K][i-1],wn);
	for(int i=K-1;i;i--)for(int j=0;j<(1<<(i-1));j++) w[i][j]=w[i+1][j<<1];
}
void NTT(poly &a, int typ=1){
	for(int i=0; i<up; i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=1,l=1; i<up; i<<=1,++l)
	for(int j=0; j<up; j+=(i<<1))
	for(int k=0; k<i; k++){
		int x=a[k+j], y=mul(w[l][k],a[k+j+i]);
		a[k+j]=add(x,y); a[k+j+i]=dec(x,y);
	} 
	if(typ==-1){
		reverse(a.begin()+1,a.end());
		for(int i=0; i<up; i++) Mul(a[i],iv[up]); 
	}
}
poly operator * (poly a, poly b){
	int deg=a.size()+b.size()-1; init(deg);
	a.resize(up); b.resize(up); NTT(a); NTT(b);
	for(int i=0; i<up; i++) Mul(a[i],b[i]); NTT(a,-1); 
	a.resize(deg); return a;
}
poly inv(poly a, int deg){
	poly b(1,ksm(a[0],Mod-2)),c;
	for(int lim=4; lim<(deg<<2); lim<<=1){
		c.resize(lim>>1); init(lim);
		for(int i=0; i<(lim>>1); i++) c[i]=i<(int)a.size()?a[i]:0;
		b.resize(up); c.resize(up); NTT(b); NTT(c);
		for(int i=0; i<up; i++) Mul(b[i],dec(2,mul(b[i],c[i])));
		NTT(b,-1); b.resize(lim>>1);
	} b.resize(deg); return b;
}
poly integ(poly a){
	a.pb(0);
	for(int i=a.size()-1;i;i--) a[i]=mul(iv[i],a[i-1]);
	a[0]=0; return a;
}
poly deriv(poly a){
	for(int i=0; i+1<(int)a.size(); i++) a[i]=mul(a[i+1],i+1);
	a.pop_back(); return a;
}
poly ln(poly a, int deg){
	a=integ(inv(a,deg)*deriv(a)); 
	a.resize(deg); return a;
}
poly Exp(poly a, int deg){
	poly b(1,1), c;
	for(int lim=2; lim<(deg<<1); lim<<=1){
		c=ln(b,lim); Dec(c[0],1);
		for(int i=0; i<lim; i++) c[i]=dec(i<(int)a.size()?a[i]:0,c[i]);
		b=b*c; b.resize(lim);
	} b.resize(deg); return b;
}
poly C, D;
void input(){
	int ma=read(), mb=read(); poly A(n+1,0), B(n+1,0); 
	while(ma--){ int x=read(); A[x]=ifac[x]; }
	while(mb--){ int x=read(); B[x]=ifac[x]; }
	B[0]=1; D.resize(n+1); C.resize(n);
	for(int i=0; i<=n; i++) D[i]=dec(B[i],A[i]); 
	for(int i=0; i<n; i++) C[i]=(i&1)?ifac[i]:Mod-ifac[i];
	C=C*A; C.resize(n+1);
}
poly Newton(int lim){
	if(lim==1) return poly(1,0);
	poly f0=Newton((lim+1)>>1);
	poly H=C; H.resize(lim); H=H*Exp(f0,lim);
	H.resize(lim); poly v=Exp(integ(H),lim);
	Dec(f0[0],1); H=H*f0; H.resize(lim);
	for(int i=0;i<lim;i++) Add(H[i],D[i]);
	poly f=integ(v*H); f.resize(lim);
	f=f*inv(v,lim); f.resize(lim);
	return f;
}
int main(){
	n=read(); 
	pre_work(n+5); NTT_init();
	input(); poly f=Newton(n+1);
	cout<<mul(fac[n],f[n]); return 0;
}