[多项式算法](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计算即可。
好像挺简单的?以前都没敢学来着
代码
这里我使用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还是比较简单易懂的。
其中还有一个小优化:对于递归到范围较小的序列,直接暴力卷积,常数更小。
参考资料: