传送门
一般当牛顿迭代里面套了Exp的时候,能写分治FFT就写分治FFT,多半跑得比牛顿迭代快。 ——zxyoi
分治FFT题解:https://blog.csdn.net/zxyoi_dreamer/article/details/101037455
题解:
首先将问题转化为如下形式:
对满足如下条件的 n n n个点的树计数,其父亲儿子的标号满足堆性质,非叶节点有 c + 2 c+2 c+2个儿子,其中 c c c个是叶子,且 c ∈ A c\in A c∈A,A集合由数据给出,另外两个儿子允许为叶子,如果不为叶子,它们的子树也要满足上述性质。(递归定义)
设 f i f_i fi表示 i i i个点的时候的答案,枚举那两个子树的大小,可以得到:
f i = 1 2 ∑ j ∑ k [ i − 1 − j − k ∈ A ] ( i − 1 j ) ( j − 1 − j k ) f j f k f_i=\frac{1}{2}\sum_{j}\sum_{k}[i-1-j-k\in A]{i-1\choose j}{j-1-j\choose k}f_jf_k fi=21j∑k∑[i−1−j−k∈A](ji−1)(kj−1−j)fjfk
注意由于两个子树是等价的,所有左右分配标号的方式我们都重复算了一遍需要除以 2 2 2,递推边界为 f 0 = 0 , f 1 = 1 f_0=0,f_1=1 f0=0,f1=1
设 f i ′ = f i i ! f_i'=\frac{f_i}{i!} fi′=i!fi, a i = [ i ∈ A ] i ! a_i=\frac{[i\in A]}{i!} ai=i![i∈A],我们可以得到一个转移: f i ′ = 1 2 i ∑ j a j ∑ k + t = i − j − 1 f k ′ f t ′ f_i'=\frac{1}{2i}\sum_ja_j\sum_{k+t=i-j-1}f_k'f_t' fi′=2i1∑jaj∑k+t=i−j−1fk′ft′。
设
f
i
f_i
fi和集合
A
A
A的EGF分别为
F
(
x
)
,
A
(
x
)
F(x),A(x)
F(x),A(x),根据上面的转移式子得到这个常微分方程:
F
′
(
x
)
=
1
2
F
2
(
x
)
A
(
x
)
+
1
F'(x)=\frac{1}{2}F^2(x)A(x)+1
F′(x)=21F2(x)A(x)+1
可以用分治NTT搞,好写还比牛顿迭代快
我们知道 A A A是一个常多项式,现在要解 F F F。
设 h ( t ) = 1 2 A ( x ) t + 1 h(t)=\frac{1}{2}A(x)t+1 h(t)=21A(x)t+1。我们要解决的问题就是 F ′ ≡ h ( F ) ( m o d x n ) F'\equiv h(F)\pmod {x^n} F′≡h(F)(modxn)
假设现在我们已经解得了 F 0 ′ ≡ h ( F 0 ) ( m o d x n 2 ) F_0'\equiv h(F_0)\pmod {x^\frac{n}{2}} F0′≡h(F0)(modx2n),考虑怎么扩展到 m o d x n \mod x^n modxn
首先将 h h h在 F 0 F_0 F0处泰勒展开,具体过程我跳了,得到:
F ′ ≡ h ( F 0 ) + h ′ ( F 0 ) F − h ′ ( F 0 ) F 0 ( m o d x n ) F'\equiv h(F_0)+h'(F_0)F-h'(F_0)F_0\pmod {x^n} F′≡h(F0)+h′(F0)F−h′(F0)F0(modxn)
但是还是没法做,看一个神奇的构造,设
v
(
x
)
=
exp
(
−
∫
h
′
(
F
0
(
x
)
)
d
x
)
v(x)=\exp(-\int h'(F_0(x)) \mathrm{d} x)
v(x)=exp(−∫h′(F0(x))dx),显然有
v
′
=
−
h
′
(
F
0
)
v
v'=-h'(F_0)v
v′=−h′(F0)v,好像是解常微分的方程可以用的套路?这是我的第一道常微分方程啊。。。
把上面的等式两边同时乘上 v v v,开始我们神奇的转化,(以下讨论全部 m o d x n \mod x^n modxn):
F ′ v = v ( h ( F 0 ) + h ′ ( F 0 ) ( F − F 0 ) ) F ′ v − h ′ ( F 0 ) v F = v ( h ( F 0 ) − h ′ ( F 0 ) F 0 ) F ′ v + v ′ F = v ( h ( F 0 ) − h ′ ( F 0 ) F 0 ) ( F v ) ′ = v ( h ( F 0 ) − h ′ ( F 0 ) F 0 ) F = 1 v ∫ v ( h ( F 0 ) − h ′ ( F 0 ) F 0 ) d x \begin{aligned} F'v&=v(h(F_0)+h'(F_0)(F-F_0))\\ F'v-h'(F_0)vF&=v(h(F_0)-h'(F_0)F_0)\\ F'v+v'F&=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{d} x \end{aligned} F′vF′v−h′(F0)vFF′v+v′F(Fv)′F=v(h(F0)+h′(F0)(F−F0))=v(h(F0)−h′(F0)F0)=v(h(F0)−h′(F0)F0)=v(h(F0)−h′(F0)F0)=v1∫v(h(F0)−h′(F0)F0)dx
发现我们并没有在意 h h h实际上是什么,也就是说这个构造套路可以用在一般的常微分方程上面。
对于这道题,我们把 h h h带入得到:
v ( x ) = exp ( − ∫ A ( x ) F 0 ( x ) d x ) F ( x ) = 1 v ( x ) ∫ v ( x ) ( 1 − 1 2 A ( x ) F 0 2 ( x ) ) d x \begin{aligned} v(x)&=\exp(-\int A(x)F_0(x)\mathrm{d}x)\\ F(x)&=\frac{1}{v(x)}\int v(x)(1-\frac{1}{2}A(x)F_0^2(x))\mathrm{d}x \end{aligned} v(x)F(x)=exp(−∫A(x)F0(x)dx)=v(x)1∫v(x)(1−21A(x)F02(x))dx
愉快迭代即可。复杂度
O
(
n
log
n
)
O(n\log n)
O(nlogn),常数根本没眼看,反正写丑了就过不了
不过vector封装多项式还挺快的,最慢的点1.6s跑出来了。
代码:
#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const
using std::cerr;
using std::cout;
cs int mod=998244353,inv2=(mod+1)>>1;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
return res;
}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}
typedef std::vector<int> Poly;
cs int bit=19,SIZE=1<<bit|1;
int r[SIZE],*w[bit+1];
int fac[SIZE],ifac[SIZE],inv[SIZE];
inline void init_NTT(){
for(int re i=1;i<=bit;++i)w[i]=new int[1<<i-1];
int wn=power(3,mod-1>>bit);w[bit][0]=1;
for(int re i=1;i<(1<<bit-1);++i)w[bit][i]=mul(w[bit][i-1],wn);
for(int re i=bit-1;i;--i)
for(int re j=0;j<(1<<i-1);++j)w[i][j]=w[i+1][j<<1];
inv[0]=inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
for(int re i=2;i<SIZE;++i){
fac[i]=mul(fac[i-1],i);
inv[i]=mul(inv[mod%i],mod-mod/i);
ifac[i]=mul(ifac[i-1],inv[i]);
}
}
inline void NTT(int *A,int len,int typ){
for(int re i=1;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
for(int re i=1,d=1;i<len;i<<=1,++d)
for(int re j=0;j<len;j+=i<<1)
for(int re k=0;k<i;++k){
int &t1=A[j+k],&t2=A[i+j+k],t=mul(t2,w[d][k]);
t2=dec(t1,t),Inc(t1,t);
}
if(typ==-1){
std::reverse(A+1,A+len);
for(int re i=0,iv=inv[len];i<len;++i)Mul(A[i],iv);
}
}
inline void NTT(Poly &A,int len,int typ){NTT(&A[0],len,typ);}
inline void init_rev(int l){
for(int re i=1;i<l;++i)r[i]=r[i>>1]>>1|((i&1)?l>>1:0);
}
inline Poly operator*(Poly a,Poly b){
int n=a.size(),m=b.size(),l=1;
while(l<n+m-1)l<<=1;
init_rev(l);
a.resize(l),NTT(a,l,1);
b.resize(l),NTT(b,l,1);
for(int re i=0;i<l;++i)a[i]=mul(a[i],b[i]);
NTT(a,l,-1);
a.resize(n+m-1);
return a;
}
inline Poly Inv(cs Poly &a,int lim){
int n=a.size();
Poly c,b(1,power(a[0],mod-2));
for(int re l=4;(l>>2)<lim;l<<=1){
init_rev(l);
c.resize(l>>1);for(int re i=0;i<(l>>1);++i)c[i]=i<n?a[i]:0;
c.resize(l);NTT(c,l,1);
b.resize(l);NTT(b,l,1);
for(int re i=0;i<l;++i)Mul(b[i],dec(2,mul(b[i],c[i])));
NTT(b,l,-1);b.resize(l>>1);
}b.resize(lim);
return b;
}
inline Poly Deriv(Poly a){
for(int re 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.push_back(0);
for(int re i=a.size()-1;i;--i)a[i]=mul(a[i-1],inv[i]);
a[0]=0;
return a;
}
inline Poly Ln(Poly a,int lim){
a=Integ(Deriv(a)*Inv(a,lim));
a.resize(lim);
return a;
}
inline Poly Exp(cs Poly &a,int lim){
int n=a.size();Poly c,b(1,1);
for(int re i=2;(i>>1)<lim;i<<=1){
c=Ln(b,i);
for(int re j=0;j<i;++j)c[j]=dec(j<n?a[j]:0,c[j]);
Inc(c[0],1);
b=b*c;b.resize(i);
}b.resize(lim);
return b;
}
cs int N=1<<18|1;
char s[N];
int n,c[N];
inline Poly Newton(int lim){
if(lim==1)return Poly(1,0);
Poly F0=Newton(lim+1>>1);
Poly H=Poly(c,c+lim)*F0;
H.resize(lim);Poly v=Integ(H);
v.resize(lim);
for(int re i=0;i<lim;++i)v[i]=dec(0,v[i]);
v=Exp(v,lim);H=H*F0;H.resize(lim);
for(int re i=0;i<lim;++i)Mul(H[i],mod-inv2);
Inc(H[0],1);Poly F=Integ(v*H);
F.resize(lim);F=F*Inv(v,lim);
F.resize(lim);return F;
}
signed main(){
#ifdef zxyoi
freopen("reaction.in","r",stdin);
#endif
init_NTT();scanf("%d%s",&n,s);
for(int re i=0;i<=n;++i)c[i]=s[i]=='1'?ifac[i]:0;
Poly F=Newton(n+1);
for(int re i=1;i<=n;++i)cout<<mul(F[i],fac[i])<<"\n";
return 0;
}