[多项式算法](Part 5)分治FFT 学习笔记

其他多项式算法传送门:

[多项式算法](Part 1)FFT 快速傅里叶变换 学习笔记

[多项式算法](Part 2)NTT 快速数论变换 学习笔记

[多项式算法](Part 3)MTT 任意模数FFT/NTT 学习笔记

[多项式算法](Part 4)FWT 快速沃尔什变换 学习笔记

分治FFT

分治FFT是一种基于CDQ分治(不是很懂和CDQ有什么关系)的算法,主要用于在\(O(nlog^2n)\)的时间复杂度内计算以下内容:

给定数列\(g_1\sim g_{n-1}\),求\(f_0\sim f_{n-1}\),其中\(f_i\)满足

\[f_0=1\\ f_i=\sum_{j=1}^jf_{i-j}\times g_j \]

分析

我们可以借鉴分治的思想,对于当前的\(f_l\sim f_r\),设\(mid=\lfloor\frac{l+r}{2}\rfloor\)

先递归计算\(f_l\sim f_{mid}\),接着考虑前半部分对后面的贡献。

设对\(f_x(mid+1\le x\le r)\)的贡献为\(w_x\),则有

\[w_x=\sum_{i=l}^{mid}f_i\times g_{x-i} \]

为了下面推导方便,把范围扩大,设\(f_i=0(mid+1\le i\le r)\)

\[w_x=\sum_{i=l}^{x-1}f_i\times g_{x-i}\\ w_x=\sum_{i=0}^{x-l-1}f_{i+l}\times g_{x-i-l} \]

此时设\(a_i=f_{i+l},b_i=g_{i+1}\),则有:

\[w_x=\sum_{i=0}^{x-l-1}a_i\times b_{x-l-1-i} \]

那么就会发现这是一个卷积的形式,使用FFT计算即可。

好像挺简单的?以前都没敢学来着

代码

例题: Luogu P4721 【模板】分治 FFT

这里我使用NTT来实现算法,当然使用FFT也是没有问题的。

加了IO优化,代码可能有点乱?

时间复杂度 \(O(nlog^2n)\)

空间复杂度 \(O(n)\)

// luogu-judger-enable-o2
#include <cmath>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
#define rint register int
typedef long long ll;

//----- IO optimize -----
#define Getchar (p1==p2&&(p2=(p1=In)+fread(In,1,1<<20,stdin),p1==p2)?EOF:*p1++)
char In[1<<20],*p1=In,*p2=In,Ch,Out[1<<21],*Outp=Out,St[15],*Tp=St;
inline int Getint(register int x=0)
{
    while(!isdigit(Ch=Getchar));
    for(;isdigit(Ch);Ch=Getchar)x=x*10+(Ch^48);
    return x;
}
inline void Putint(int x,char c)
{
    do *Tp++=x%10^48;while(x/=10);
    do *Outp++=*--Tp;while(Tp!=St);
    *Outp++=c;
}

const int p=998244353,g1=3,g2=(p+1)/3;//g2为1/3 mod p
inline int Add(int a,int b){return (a+=b)>=p?a-p:a;}

inline ll Pow(ll a,ll b)
{
    ll Res=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)Res=Res*a%p;
    return Res%p;
}

namespace Poly
{
    int r[1<<18];
    void NTT(int n,int* A,const int g)//没学过NTT?大丈夫,换成你喜欢的FFT即可
    {
        for(rint i=0;i<n;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
        for(rint i=2,h=1;i<=n;i<<=1,h<<=1)
            for(rint j=0,Rs=Pow(g,(p-1)/i);j<n;j+=i)
                for(rint k=0,Root=1;k<h;++k,Root=(ll)Root*Rs%p)
                {
                    int Tmp=(ll)A[j+h+k]*Root%p;
                    A[j+h+k]=Add(A[j+k],p-Tmp),A[j+k]=Add(A[j+k],Tmp);
                }
    }

    int A[1<<18],B[1<<18];

    void Multiply(int n,int *A,int *B)//A=A*B
    {
        for(rint i=1,l=(int)log2(n);i<n;++i)
            r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        NTT(n,A,g1),NTT(n,B,g1);
        for(rint i=0;i<n;++i)A[i]=(ll)A[i]*B[i]%p;
        NTT(n,A,g2);
        int Invn=Pow(n,p-2);
        for(rint i=0;i<n;++i)A[i]=(ll)A[i]*Invn%p;
    }

    void CDQ_FFT(int* f,int* g,int l,int r)
    {
        if(l==r)return;
        int Mid=(l+r)>>1,Len=(r-l+1),n=1;
        CDQ_FFT(f,g,l,Mid);
        while(n<=Len)n<<=1;
        memset(A,0,n*sizeof(int));
        memset(B,0,n*sizeof(int));
        for(rint i=l;i<=Mid;++i)A[i-l]=f[i];
        for(rint i=1;i<=r-l;++i)B[i-1]=g[i];
        Multiply(n,A,B);
        for(rint i=Mid+1;i<=r;++i)f[i]=Add(f[i],A[i-l-1]);
        CDQ_FFT(f,g,Mid+1,r);
    }
}

int n,f[1<<18],g[1<<18];

int main()
{
    n=Getint(),f[0]=1;
    for(rint i=1;i<n;++i)g[i]=Getint();
    Poly::CDQ_FFT(f,g,0,n-1);
    for(rint i=0;i<n;++i)Putint(f[i],i==n-1?'\n':' ');
    return fwrite(Out,1,Outp-Out,stdout),0;
}

分治FFT还是比较简单易懂的。

其中还有一个小优化:对于递归到范围较小的序列,直接暴力卷积,常数更小。

参考资料:

分治FFT 学习笔记-igronemyk

posted @ 2019-08-05 20:51  LanrTabe  阅读(2015)  评论(0编辑  收藏  举报