分治FFT

更新日志

2025/01/09：开工。

2025/01/09：使用封装模板重写代码。

---

概念

借助分治的思想，在 \(O(n\log^2n)\) 的复杂度内计算下述问题：

给出 \(G(x),x\in[1,n)\)，求 \(F(x),x\in[0,n)\)，\(F(x)\) 满足：

\[F(0)=0\\ \forall x\in(0,n),F(x)=\sum_{i=1}^xF(x-i)G(i) \]

思路

我们借鉴分治思想，若当前操作区间为 \([l,r]\)，令 \(mid=\lfloor\frac{l+r}{2}\rfloor\)。

我们先计算出左区间内部的贡献，现在考虑整个左区间对右区间的贡献。

若对 \(F(x)\) 的贡献为 \(w_x\)（\(x\in(mid,r]\)）：

\[w_x=\sum_{i=l}^{mid}F(i)G(x-i) \]

不妨令此时 \(\forall x\in(mid,r],F(x)=0\)，那么这个式子可以改写成：

\[\begin{align*} w_x=&\sum_{i=l}^{x-1}F(i)G(x-i)\\ =&\sum_{i=0}^{x-l-1}F(i+l)G(x-l-i) \end{align*} \]

设 \(t=x-l-1,A(x)=F(x+l),B(x)=G(x+1)\)，则：

\[w_x=\sum_{i=0}^tA(i)B(t-i) \]

这就是一个卷积，\(\text{FFT}/\text{NTT}\) 求解即可。

例题代码

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double db;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define repl(i,x,y) for(int i=(x);i<(y);i++)
#define per(i,x,y) for(int i=(x);i>=(y);i--)

const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;

namespace NTT{
    const int p=998244353,g=3,gi=332748118;
    typedef vec<ll> poly;
    ll qpow(ll a,ll b){
        ll res=1;
        while(b){
            if(b&1)(res*=a)%=p;
            (a*=a)%=p;
            b>>=1;
        }
        return res;
    }

    vec<int> r;
    void ntt(int lim,poly &a,int type){
        a.resize(lim);
        repl(i,0,lim)if(i<r[i])swap(a[i],a[r[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            ll w1=qpow((~type)?g:gi,(p-1)/(mid<<1));
            for(int j=0;j<lim;j+=(mid<<1)){
                ll w=1;
                repl(k,0,mid){
                    ll x=a[j+k],y=w*a[j+mid+k]%p;
                    a[j+k]=(x+y)%p;
                    a[j+mid+k]=(x-y+p)%p;
                    w=w*w1%p;
                }
            }
        }
    }
    poly operator*(poly a,poly b){
        int lim=1,l=0,len=a.size()+b.size()-1;
        while(lim<=a.size()+b.size())lim<<=1,l++;
        r.resize(lim);
        repl(i,0,lim)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        ntt(lim,a,1);ntt(lim,b,1);
        repl(i,0,lim)a[i]=(a[i]*b[i])%p;
        ntt(lim,a,-1);
        a.resize(len);
        ll inv=qpow(lim,p-2);
        repl(i,0,a.size())(a[i]*=inv)%=p;
        return a;
    }
    poly operator+(poly a,poly b){
        if(a.size()<b.size())a.resize(b.size());
        repl(i,0,b.size())a[i]=(a[i]+b[i])%p;
        return a;
    }
    poly operator-(poly a,poly b){
        if(a.size()<b.size())a.resize(b.size());
        repl(i,0,b.size())a[i]=(a[i]-b[i]+p)%p;
        return a;
    }
}
using namespace NTT;

const int N=1e5+5;

int n;
ll gg[N];
ll f[N];

void solve(int l,int r){
    if(l==r)return;
    int m=l+r>>1;
    solve(l,m);
    poly a(r-l),b(r-l);
    rep(i,l,m)a[i-l]=f[i];
    repl(i,0,r-l)b[i]=gg[i+1];
    poly c=a*b;
    rep(x,m+1,r)(f[x]+=c[x-l-1])%=p;
    solve(m+1,r);
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin>>n;
    repl(i,1,n)cin>>gg[i];
    f[0]=1;
    solve(0,n-1);
    repl(i,0,n)cout<<f[i]<<" ";
    return 0;
}