题解 P4726 【【模板】多项式指数函数】

cosmicAC · 2019-01-16 20:03:34 · 题解

发一个分治NTT的题解。如果想要看一个log的正解请参考其他dalao的代码。

题目要求的是B(x)=e^{A(x)} 两边同时求导得到B'(x)=A'(x)e^{A(x)} 即B'(x)=A'(x)B(x) 同时积分得到B(x)=\int{A'(x)B(x)} 又因为将x=0代入A和B中可以得到B(0)=1，可以以此为初始值进行分治FFT，每次把A'和B的卷积平移一项后加到B自身上。做多项式积分时第n项会产生\frac{1}{n}的常数，可以在递归边界时乘上去。

时间复杂度O(n\log^2{n})，但因为常数较小，实际上只比牛顿迭代慢一倍以内。

代码如下（我的NTT写法比较短）：

// luogu-judger-enable-o2
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll lg=18,N=1<<lg,mod=998244353;
int n;
ll f[N],g[N],rev[N],w[N],c[N],d[N];
ll po(ll a,ll b){
    ll r=1;for(;b;b>>=1,(a*=a)%=mod)if(b&1)(r*=a)%=mod;
    return r;
}
void FFT(ll *a,int lg){
    for(int i=0;i<(1<<lg);i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=0;i<lg;i++)
        for(int j=0;j<(1<<lg);j++)if(j&1<<i){
            ll x=j^1<<i,l=a[x],r=a[j]*w[(j&(1<<i)-1)<<lg-1-i]%mod;
            a[x]=(l+r)%mod,a[j]=(l+mod-r)%mod;
        }
}
void work(int l,int r){ //此处是左闭右开区间[l,r)
    if(l+1==r){
        if(l)(f[l]*=po(l,mod-2))%=mod;else f[l]=1;
        return;
    }
    int mid=l+r>>1,_lg,n;
    work(l,mid);
    for(_lg=0;1<<_lg<r-l-1;_lg++);
    n=1<<_lg;
    for(int i=0;i<n;i++)
        c[i]=d[i]=0,w[i]=po(3,1ll*(mod-1)*i/n),
        rev[i]=rev[i/2]/2|(i&1)<<_lg-1;
    for(int i=0;i<mid-l;i++)c[i]=f[i+l];
    for(int i=0;i<r-l-1;i++)d[i]=g[i];
    FFT(c,_lg),FFT(d,_lg);
    for(int i=0;i<n;i++)(c[i]*=d[i])%=mod,w[i]=po(w[i],mod-2);
    FFT(c,_lg);
    for(int i=mid-1-l;i<r-l-1;i++)(f[i+l+1]+=c[i]*po(n,mod-2)%mod)%=mod;
    work(mid,r);
}
int main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%lld",g+i);
    for(int i=1;i<n;i++)g[i-1]=g[i]*i%mod;g[n-1]=0;
    work(0,n);
    for(int i=0;i<n;i++)printf("%lld%c",f[i],i==n-1?'\n':' ');
    return 0;
}