可爱的分块遇见可爱的 FFT

· · 题解

由于 \text{FFT} 太可爱了，交一发直接最优解 \text{(2022.9.30)}，比目前所有 \text{NTT} 都要快。

对于所有模数为 998244353 的 \text{NTT} 都可以被如下数据卡掉（仅限模数为 998244353）：

100000
1000 1000/*100000个“1000”*/
1000000
1 100000 0
1 100000 0
/*100000行“1 100000 0”*/

这样单次 \text{NTT} 的结果会达到 10^9，当然一般不会有人去卡，这样没意义。

这题一看就是一个加法卷积形式，无奈不可能每一次操作都卷一下，这样是 O(nm\log_2n) 的。

并没有一种方法可以将多个区间卷单点一起做，但如果区间一样就不同了，因为对于完全相同的区间，我们完全可以将操作加到一个多项式，在和该区间做卷积。

我们可以进行分块，对于单块的操作一起算，设有 K 个块，则复杂度为 O(K(n\log_2n+m)+\frac{nm}{K}，令两式相等，解得 K=\sqrt{\frac{nm}{n\log_2n+m}}，此时时间复杂度为 \sqrt{nm(n\log_2n+m)}，取 m=n\log_2n，则时间复杂度为 n\sqrt n\log_2n，在 n=10^5 时有点紧。

具体的，处理区间 [bl,br] 时，可以将 b_{L+k-l}+=a_k 转化为 b_{L+k-l}=a_kf_{L-l}，而 f 可以扫描所有整块询问累加，考虑到可能的 L<l 的情况，可以将 k 右移 bl，L-l 左移 bl，保证下标均在 [1,n] 之间。

十分开心地做完了：

#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5,M=1e6+6;
const double pai=3.1415926535897932;
int a[N],b[N],n,q,m;
int l_[M],L[M],R[M],tr[M];
namespace fast_io{
    char bf[N+5],ob[N+38];
    int it,ed,f,c,ot,t,stk[38],x;
#define gc (it==ed&&(ed=(it=0)+fread(bf,1,N,stdin),it==ed)?EOF:bf[it++])
    inline int read(){
        for(c=gc;c<48;c=gc);
        for(x=0;c>47;x=x*10+(48^c),c=gc);
        return x;
    }inline void fls(){
        fwrite(ob,1,ot,stdout),ot=0;
    }inline void write(int x){
        for(t=0;x>9;stk[++t]=48^(x%10),x/=10);
        for(ob[ot++]=48^x;t;ob[ot++]=stk[t--]);
        ob[ot++]='\n';if(ot>N)fls();
    }
}using fast_io::read;
using fast_io::write;
struct cp{
    double x,y;
    inline void operator*=(cp &z){
        *this={x*z.x-y*z.y,x*z.y+y*z.x};
    }inline cp operator*(cp &z){
        return {x*z.x-y*z.y,x*z.y+y*z.x};
    }inline void operator/=(int &p){
        x/=p,y/=p;
    }inline cp operator+(cp &z){
        return {x+z.x,y+z.y};
    }inline cp operator-(cp &z){
        return {x-z.x,y-z.y};
    }inline cp operator+=(int p){
        x+=p;
    }
}f[M],nw,rt,fl,fr;
inline void fft(bool tp){
    if(tp)reverse(f+1,f+m);
    int ln,md,r,i;
    for(i=0;i<m;++i)
        if(i<tr[i])swap(f[i],f[tr[i]]);
    for(md=1;md<m;md<<=1){
        r=md<<1,rt={cos(pai/md),sin(pai/md)};
        for(ln=0;ln<m;ln+=r)
            for(i=0,nw={1};i<md;++i,nw*=rt){
                fl=f[ln|i],fr=nw*f[ln|md|i];
                f[ln|i]=fl+fr,f[ln|md|i]=fl-fr;
            }
    }if(tp)for(i=0;i<m;++i)f[i]/=m;
}
int main(){
    int i,l,r,_l,_r,k;n=read();
    for(l=1;l<=n;++l)a[l]=read();q=read();
    for(i=1;i<=q;++i)L[i]=read(),R[i]=read(),l_[i]=read();
    for(m=2,k=0;m<=n;m<<=1,++k);
    for(i=1;i<m;++i)tr[i]=(tr[i>>1]>>1)|((i&1)<<k);
    for(l=1;l<=n;l=r+1){
        r=min(l+2048,n);
        for(i=0;i<m;++i)f[i]={0,0};
        for(i=1;i<=q;++i)
            if(L[i]<=l&&R[i]>=r)f[l_[i]+l-L[i]].y+=1;
            else if(L[i]<=r&&R[i]>=l){
                _l=max(l,L[i]),_r=min(r,R[i]);
                for(k=_l;k<=_r;++k)b[l_[i]+k-L[i]]+=a[k];
            }
        for(i=l;i<=r;++i)f[i-l].x+=a[i];
        for(i=0,fft(0);i<m;++i)f[i]*=f[i];
        for(i=1,fft(1);i<=n;++i)b[i]+=(f[i].y+1)/2;
    }for(i=1;i<=n;++i)write(b[i]);
    fast_io::fls();return 0;
}