UOJ#50 【UR#3C】链式反应 FFT求解多项式线性常微分方程

链式反应 FFT求解多项式线性常微分方程"/>

UOJ#50 【UR#3C】链式反应 FFT求解多项式线性常微分方程

题目大意：给定 n $n$和集合C$C$，对于 i=1..n $i=1..n$求多少 i $i$个节点有标号的多叉树满足：
1.父亲节点的标号大于子节点
2.一个点如果有儿子，则有两个无序的α$\alpha$型儿子，有 c $c$个无序的β$\beta$型儿子，其中 c∈C $c\in C$
3.如果一个点是根节点或 α $\alpha$型儿子，那么它可以有儿子或者是一个叶节点；如果一个点是 β $\beta$型儿子，那么它只能是一个叶节点

由于有标号，所以这里显然要使用指数级生成函数
设 fi $f_i$表示 i $i$个节点的多叉数数量，其指数级生成函数为：
F(x)=∑fixii!$F(x)=\sum f_i\frac{x^i}{i!}$
那么考虑：
根节点是固定的，不参与标号的排列，首先把根节点刨掉，即：
F′(x)=∑fixi−1(i−1)! $F'(x)=\sum f_i\frac{x^{i-1}}{(i-1)!}$
然后一棵树可能只有一个叶节点，刨掉之后就什么都不剩了；
也可能有儿子，这时候儿子分为三部分：
α1,α2,β $\alpha1,\alpha2,\beta$，三者的指数级生成函数分别为 F(x),F(x),C(x) $F(x),F(x),C(x)$
故直观上来看应该是 F2(x)C(x) $F^2(x)C(x)$，但是这样不对，因为 α1 $\alpha1$和 α2 $\alpha2$没有顺序，所以应该是 12F2(x)C(x) $\frac12F^2(x)C(x)$
可以列出方程：
F′(x)=12F2(x)C(x)+1 $F'(x)=\frac12F^2(x)C(x)+1$

这个方程怎么解？

牛顿迭代。

假设现在我们有待定多项式 x(t) $x(t)$、常多项式 a(t) $a(t)$、已知函数 f(x)=12ax2+1 $f(x)=\frac12ax^2+1$，求解方程：
ddtx=f(x) $\frac d{dt}x=f(x)$
老办法，倍增处理。
假设我们已经知道了 x(t) $x(t)$的前 n $n$项xn$x_n$，求 x(t) $x(t)$的前 2n $2n$项 x2n $x_{2n}$
泰勒展开：
ddtx2n=f(x2n)=f(xn)+f′(xn)(x2n−xn)+f′′(xn)(x2x−xn)2+... $\frac d{dt}x_{2n}=f(x_{2n})=f(x_n)+f'(x_n)(x_{2n}-x_n)+f''(x_n)(x_{2x}-x_n)^2+...$
ddtx2n=f(xn)+f′(xn)(x2n−xn) (mod t2n) $\frac d{dt}x_{2n}=f(x_n)+f'(x_n)(x_{2n}-x_n)\ (mod\ t^{2n})$
ddtx2n−x2nf′(xn)=f(xn)−xnf′(xn) (mod t2n) $\frac d{dt}x_{2n}-x_{2n}f'(x_n)=f(x_n)-x_nf'(x_n)\ (mod\ t^{2n})$
看到这里学过微积分的人已经会构造了吧。
令 r=e−∫f′(xn)dt $r=e^{-\int f'(x_n)dt}$，则 ddtr=−f′(xn)r $\frac d{dt}r=-f'(x_n)r$
等式两边同乘 r $r$得到：
rddtx2n+x2nddtr=(f(xn)−xnf′(xn))r (mod t2n)$r\frac d{dt}x_{2n}+x_{2n}\frac d{dt}r=(f(x_n)-x_nf'(x_n))r\ (mod\ t^{2n})$
ddt(x2nr)=(f(xn)−xnf′(xn))r (mod t2n) $\frac d{dt}(x_{2n}r)=(f(x_n)-x_nf'(x_n))r\ (mod\ t^{2n})$
x2n=∫(f(xn)−xnf′(xn))r dtr (mod t2n) $x_{2n}=\frac{\int(f(x_n)-x_nf'(x_n))r\ dt}r\ (mod\ t^{2n})$
x2n=∫(1−12ax2n)r dtr (mod t2n) $x_{2n}=\frac{\int(1-\frac12ax_n^2)r\ dt}r\ (mod\ t^{2n})$
里面的一切计算都是 O(nlogn) $O(nlogn)$的，总时间复杂度 T(n)=O(nlogn)+T(n2)=O(nlogn) $T(n)=O(nlogn)+T(\frac n2)=O(nlogn)$
注意常数……

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 530000
#define MOD 998244353
#define G 3
using namespace std;
int n,m,d;
long long inv[M];
int A[M],B[M];
void Linear_Shaker()
{int i;for(inv[1]=1,i=2;i<=d<<1;i++)inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
}
long long Quick_Power(long long x,int y)
{long long re=1;while(y){if(y&1) (re*=x)%=MOD;(x*=x)%=MOD; y>>=1;}return re;
}
void FFT(int a[],int n,int type){static int rev_bit[M];int i,j,k,w,wn,t,bit;for(bit=0;1<<bit<n;++bit);static int b[M];for(i=0;i<n;i++)rev_bit[i]=(rev_bit[i>>1]>>1)|((i&1)<<bit-1);for(i=0;i<n;i++)if(i<rev_bit[i])swap(a[i],a[rev_bit[i]]);for(k=2;k<=n;k<<=1)for(wn=Quick_Power(G,(long long)(MOD-1)/k*type%(MOD-1)),i=0;i<n;i+=k)for(w=1,j=0;j<k>>1;++j,w=(long long)w*wn%MOD){t=(long long)w*a[i+j+(k>>1)]%MOD;a[i+j+(k>>1)]=a[i+j]-t<0?a[i+j]-t+MOD:a[i+j]-t;a[i+j]=a[i+j]+t>=MOD?a[i+j]+t-MOD:a[i+j]+t;}if(type!=1){for(i=0;i<n;i++)a[i]=(long long)a[i]*inv[n]%MOD;}
}
void Get_Inv(int a[],int b[],int n) 
//求a的逆，长度为n，结果储存在b中。
//要求传入时b[0...n<<1]为空。
{static int temp[M];int i;if(n==1){b[0]=Quick_Power(a[0],MOD-2);return ;}Get_Inv(a,b,n>>1);memcpy(temp,a,sizeof(a[0])*n);memset(temp+n,0,sizeof(a[0])*n);FFT(temp,n<<1,1);FFT(b,n<<1,1);for(i=0;i<n<<1;i++)temp[i]=(long long)b[i]*(2-(long long)temp[i]*b[i]%MOD+MOD)%MOD;FFT(temp,n<<1,MOD-2);memcpy(b,temp,sizeof(a[0])*n);memset(b+n,0,sizeof(a[0])*n);
}
void Get_Ln(int a[],int b[],int n)
//求a的ln，长度为n，结果储存在b中。
//要求a的常数项为1。
{static int a_[M],a_inv[M];int i;Get_Inv(a,a_inv,n);for(i=0;i<n-1;i++)a_[i]=(long long)a[i+1]*(i+1)%MOD;FFT(a_,n<<1,1);FFT(a_inv,n<<1,1);for(i=0;i<n<<1;i++)b[i]=(long long)a_[i]*a_inv[i]%MOD;FFT(b,n<<1,MOD-2);for(i=n-1;i;i--)b[i]=b[i-1]*inv[i]%MOD;b[0]=0;memset(b+n,0,sizeof(b[0])*n);memset(a_,0,sizeof(a_[0])*n<<1);memset(a_inv,0,sizeof(a_inv[0])*n<<1);
}
void Get_Exp(int a[],int b[],int n)
//求a的exp，长度为n，结果储存在b中。
//要求传入时b[0...n<<1]为空。
//要求a的常数项为0。
{static int temp[M];int i;if(n==1){b[0]=1;return ;}Get_Exp(a,b,n>>1);Get_Ln(b,temp,n);for(i=0;i<n;i++)temp[i]=((i==0)+MOD-temp[i]+a[i])%MOD;FFT(temp,n<<1,1);FFT(b,n<<1,1);for(i=0;i<n<<1;i++)b[i]=(long long)b[i]*temp[i]%MOD;FFT(b,n<<1,MOD-2);memset(b+n,0,sizeof(b[0])*n);
}
void Newton_Method(int a[],int b[],int n)
//求解常微分方程dx/dt=ax^2/2+1，长度为n，结果储存在b中。
//要求传入时b[0...n<<1]为空。
//x[2n]=int( (1-ax^2/2)*r )/r
//r=exp(-int(ax))
{static int temp[M],temp_temp[M],r[M],r_temp[M];int i;if(n==1){b[0]=0;return ;}Newton_Method(a,b,n>>1);memcpy(temp,a,sizeof(a[0])*n);memset(temp+n,0,sizeof(a[0])*n);FFT(temp,n<<1,1);FFT(b,n<<1,1);for(i=0;i<n<<1;i++)temp_temp[i]=(long long)temp[i]*b[i]%MOD;FFT(temp_temp,n<<1,MOD-2);for(i=n-1;i;i--)temp_temp[i]=temp_temp[i-1]*inv[i]%MOD*(MOD-1)%MOD;temp_temp[0]=0;memset(r,0,sizeof(a[0])*n<<1);Get_Exp(temp_temp,r,n);for(int i=0;i<n<<1;i++)temp[i]=(1+inv[2]*(MOD-1)%MOD*temp[i]%MOD*b[i]%MOD*b[i]%MOD)%MOD;FFT(temp,n<<1,MOD-2);memset(temp+n,0,sizeof(a[0])*n);memcpy(r_temp,r,sizeof(a[0])*n<<1);FFT(temp,n<<1,1);FFT(r_temp,n<<1,1);for(int i=0;i<n<<1;i++)temp[i]=(long long)temp[i]*r_temp[i]%MOD;FFT(temp,n<<1,MOD-2);for(i=n-1;i;i--)temp[i]=temp[i-1]*inv[i]%MOD;temp[0]=0;memset(temp+n,0,sizeof(a[0])*n);memset(r_temp,0,sizeof(a[0])*n<<1);Get_Inv(r,r_temp,n);FFT(temp,n<<1,1);FFT(r_temp,n<<1,1);for(int i=0;i<n<<1;i++)b[i]=(long long)temp[i]*r_temp[i]%MOD;FFT(b,n<<1,MOD-2);memset(b+n,0,sizeof(a[0])*n);memset(temp+n,0,sizeof(a[0])*n);
}
int main()
{cin>>n;for(d=1;d<=n+1;d<<=1);Linear_Shaker();long long temp=1;for(int i=0;i<n;i++){scanf("%1d",&A[i]);A[i]=A[i]*temp;(temp*=inv[i+1])%=MOD;}Newton_Method(A,B,d);temp=1;for(int i=1;i<=n;i++){(temp*=i)%=MOD;printf("%d\n",int(B[i]*temp%MOD));}return 0;
}