UOJ#50 【UR#3C】鍊式反應 FFT求解多項式線性常微分方程

題目大意：給定 n 和集合C，對于 i=1..n 求多少 i 個節點有标号的多叉樹滿足：

1.父親節點的标号大于子節點

2.一個點如果有兒子，則有兩個無序的α型兒子，有 c 個無序的β型兒子，其中 c∈C

3.如果一個點是根節點或 α 型兒子，那麼它可以有兒子或者是一個葉節點；如果一個點是 β 型兒子，那麼它隻能是一個葉節點

由于有标号，是以這裡顯然要使用指數級生成函數

設 fi 表示 i 個節點的多叉數數量，其指數級生成函數為：

F(x)=∑fixii!

那麼考慮：

根節點是固定的，不參與标号的排列，首先把根節點刨掉，即：

F′(x)=∑fixi−1(i−1)!

然後一棵樹可能隻有一個葉節點，刨掉之後就什麼都不剩了；

也可能有兒子，這時候兒子分為三部分：

α1,α2,β ，三者的指數級生成函數分别為 F(x),F(x),C(x)

故直覺上來看應該是 F2(x)C(x) ，但是這樣不對，因為 α1 和 α2 沒有順序，是以應該是 12F2(x)C(x)

可以列出方程：

F′(x)=12F2(x)C(x)+1

這個方程怎麼解？

牛頓疊代。

假設現在我們有待定多項式 x(t) 、常多項式 a(t) 、已知函數 f(x)=12ax2+1 ，求解方程：

ddtx=f(x)

老辦法，倍增處理。

假設我們已經知道了 x(t) 的前 n 項xn，求 x(t) 的前 2n 項 x2n

泰勒展開：

ddtx2n=f(x2n)=f(xn)+f′(xn)(x2n−xn)+f′′(xn)(x2x−xn)2+...

ddtx2n=f(xn)+f′(xn)(x2n−xn) (mod t2n)

ddtx2n−x2nf′(xn)=f(xn)−xnf′(xn) (mod t2n)

看到這裡學過微積分的人已經會構造了吧。

令 r=e−∫f′(xn)dt ，則 ddtr=−f′(xn)r

等式兩邊同乘 r 得到：

rddtx2n+x2nddtr=(f(xn)−xnf′(xn))r (mod t2n)

ddt(x2nr)=(f(xn)−xnf′(xn))r (mod t2n)

x2n=∫(f(xn)−xnf′(xn))r dtr (mod t2n)

x2n=∫(1−12ax2n)r dtr (mod t2n)

裡面的一切計算都是 O(nlogn) 的，總時間複雜度 T(n)=O(nlogn)+T(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[]=,i=;i<=d<<;i++)
        inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
}
long long Quick_Power(long long x,int y)
{
    long long re=;
    while(y)
    {
        if(y&) (re*=x)%=MOD;
        (x*=x)%=MOD; y>>=;
    }
    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=;<<bit<n;++bit);
    static int b[M];
    for(i=;i<n;i++)
        rev_bit[i]=(rev_bit[i>>]>>)|((i&)<<bit-);
    for(i=;i<n;i++)
        if(i<rev_bit[i])
            swap(a[i],a[rev_bit[i]]);
    for(k=;k<=n;k<<=)
        for(wn=Quick_Power(G,(long long)(MOD-)/k*type%(MOD-)),i=;i<n;i+=k)
            for(w=,j=;j<k>>;++j,w=(long long)w*wn%MOD)
            {
                t=(long long)w*a[i+j+(k>>)]%MOD;
                a[i+j+(k>>)]=a[i+j]-t<?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!=)
    {
        for(i=;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==)
    {
        b[]=Quick_Power(a[],MOD-);
        return ;
    }
    Get_Inv(a,b,n>>);
    memcpy(temp,a,sizeof(a[])*n);
    memset(temp+n,,sizeof(a[])*n);
    FFT(temp,n<<,);
    FFT(b,n<<,);
    for(i=;i<n<<;i++)
        temp[i]=(long long)b[i]*(-(long long)temp[i]*b[i]%MOD+MOD)%MOD;
    FFT(temp,n<<,MOD-);
    memcpy(b,temp,sizeof(a[])*n);
    memset(b+n,,sizeof(a[])*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=;i<n-;i++)
        a_[i]=(long long)a[i+]*(i+)%MOD;
    FFT(a_,n<<,);
    FFT(a_inv,n<<,);
    for(i=;i<n<<;i++)
        b[i]=(long long)a_[i]*a_inv[i]%MOD;
    FFT(b,n<<,MOD-);
    for(i=n-;i;i--)
        b[i]=b[i-]*inv[i]%MOD;
    b[]=;
    memset(b+n,,sizeof(b[])*n);
    memset(a_,,sizeof(a_[])*n<<);
    memset(a_inv,,sizeof(a_inv[])*n<<);
}
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==)
    {
        b[]=;
        return ;
    }
    Get_Exp(a,b,n>>);
    Get_Ln(b,temp,n);
    for(i=;i<n;i++)
        temp[i]=((i==)+MOD-temp[i]+a[i])%MOD;
    FFT(temp,n<<,);
    FFT(b,n<<,);
    for(i=;i<n<<;i++)
        b[i]=(long long)b[i]*temp[i]%MOD;
    FFT(b,n<<,MOD-);
    memset(b+n,,sizeof(b[])*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==)
    {
        b[]=;
        return ;
    }
    Newton_Method(a,b,n>>);
    memcpy(temp,a,sizeof(a[])*n);
    memset(temp+n,,sizeof(a[])*n);
    FFT(temp,n<<,);
    FFT(b,n<<,);
    for(i=;i<n<<;i++)
        temp_temp[i]=(long long)temp[i]*b[i]%MOD;
    FFT(temp_temp,n<<,MOD-);
    for(i=n-;i;i--)
        temp_temp[i]=temp_temp[i-]*inv[i]%MOD*(MOD-)%MOD;
    temp_temp[]=;
    memset(r,,sizeof(a[])*n<<);
    Get_Exp(temp_temp,r,n);

    for(int i=;i<n<<;i++)
        temp[i]=(+inv[]*(MOD-)%MOD*temp[i]%MOD*b[i]%MOD*b[i]%MOD)%MOD;
    FFT(temp,n<<,MOD-);
    memset(temp+n,,sizeof(a[])*n);

    memcpy(r_temp,r,sizeof(a[])*n<<);
    FFT(temp,n<<,);
    FFT(r_temp,n<<,);
    for(int i=;i<n<<;i++)
        temp[i]=(long long)temp[i]*r_temp[i]%MOD;
    FFT(temp,n<<,MOD-);

    for(i=n-;i;i--)
        temp[i]=temp[i-]*inv[i]%MOD;
    temp[]=;
    memset(temp+n,,sizeof(a[])*n);

    memset(r_temp,,sizeof(a[])*n<<);
    Get_Inv(r,r_temp,n);

    FFT(temp,n<<,);
    FFT(r_temp,n<<,);
    for(int i=;i<n<<;i++)
        b[i]=(long long)temp[i]*r_temp[i]%MOD;
    FFT(b,n<<,MOD-);

    memset(b+n,,sizeof(a[])*n);
    memset(temp+n,,sizeof(a[])*n);
}
int main()
{
    cin>>n;
    for(d=;d<=n+;d<<=);
    Linear_Shaker();

    long long temp=;
    for(int i=;i<n;i++)
    {
        scanf("%1d",&A[i]);
        A[i]=A[i]*temp;
        (temp*=inv[i+])%=MOD;
    }

    Newton_Method(A,B,d);

    temp=;
    for(int i=;i<=n;i++)
    {
        (temp*=i)%=MOD;
        printf("%d\n",int(B[i]*temp%MOD));
    }

    return ;
}