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

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

由于有标号,所以这里显然要使用指数级生成函数
fi 表示 i 个节点的多叉数数量,其指数级生成函数为:
F(x)=fixii!
那么考虑:
根节点是固定的,不参与标号的排列,首先把根节点刨掉,即:
F(x)=fixi1(i1)!
然后一棵树可能只有一个叶节点,刨掉之后就什么都不剩了;
也可能有儿子,这时候儿子分为三部分:
α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)(x2nxn)+f′′(xn)(x2xxn)2+...
ddtx2n=f(xn)+f(xn)(x2nxn) (mod t2n)
ddtx2nx2nf(xn)=f(xn)xnf(xn) (mod t2n)
看到这里学过微积分的人已经会构造了吧。
r=ef(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=(112ax2n)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[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;
}
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值