分治FFT学习笔记

分治$\rm FFT$与其优化

---

前置知识 ：

快速傅里叶变换，快速数论变换，多项式求逆，CDQ分治，生成函数。
【（蒟蒻的学习笔记）$\rm FFT$ $\rm NTT$】【或者可以自行去看其他更好的博客讲解】

---

进入正题

• 引入问题

已知函数$g(x)$在$x\in [1,n-1]$的函数值，求给定的$f(x)$的函数在$x\in [1,n-1]$的值，其中$f(x)$定义如下：
$f(x)=\sum\limits_{y=1}^xf(x-y)g(y)$
其中，$f(0)=1$。（答案对$998244353$取模）

数据范围：$2\leq n\leq10^5$，其余均小于模数。

---

• 我会暴力！

其实，根据定义式，直接就可以写出$O(n^2)$（实际是$\frac{n\times(n-1)}{2}$）的复杂度的算法。

但是显然不行。

• 我会$\rm FFT$！

这个其实很容易看出是一个卷积式子，但是对于每一个$f(x)$，我们不可能重新的去算一次$\rm FFT$，因为那样的话复杂度会退化成$O(n^2logn)$，还不如暴力了，所以我们要考虑优化。

如果我们已经知道了$f(1)\sim f(\frac{n}{2})$，那么我们可以先计算出这些已知的对于未知部分$f(\frac{n}{2}+1)\sim f(n)$的贡献，我们看原式，发现对于一个已知的$f(x)$，对于后面的$f(i)$都只会有$f(i-j)g(i)[i-j=x]$的贡献。

所以对于一个当前的区间$l\sim r,mid=\left\lfloor\frac{l+r}{2}\right\rfloor$我们已知$f(l)\sim f(mid)$，那么它对于$f(mid+1)\sim f(r)$的贡献，可以用如下方法快速计算$O((r-l+1)log(r-l+1))$:

我们令$A(i)=f(i+l)\{i\in[0,mid-l]\}$，再令$B(i)=g(i)\{i\in[1,r-l]\}$，那么我们再令$C(i)=A(i)\bigotimes B(i)$（卷积，相当于$C(x)=\sum A(i)\times B(j)[i+j=x]$）,那么对于$f(mid+1)\sim f(r)$的其中的$f(x)$贡献就为$C(x-l-1)$，这个就是用$\rm FFT$实现了。

那么对于整个区间$[1,n]$，我们分治下去，先计算$[1,\frac{n+1}{2}]$，这时你已经知道前半部分的值，然后$\rm FFT$计算贡献，加到后半部分，然后去算$[\frac{n+1}{2},n]$的值，这也就是$\rm CDQ$分治的过程。

总的复杂度为$O(nlog^2n)$（总共递归$logn$层，每层元素$n$个，$\rm FFT$为$O(nlogn)$，所以总复杂度为$O(\text{分治}\times \text{FFT})$）。

Luogu 模板
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=4e5+10;
const ll mod=998244353ll,G=3;//模数与原根
ll inv[M],g[M],f[M];
int R[M],lg,sze;

ll fpow(ll a,ll b){
	ll ans=1;
	for(;b;b>>=1,a=(a*a)%mod){
		if(b&1) ans=(ans*a)%mod;
	}
	return ans;
}

void NTT(ll *a,int n,int f){
	for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);
	for(int i=2;i<=n;i<<=1){
		int now=i>>1;
		ll wn=inv[i];//现场算或者预处理
//		fpow(G,(mod-1)/i);
		for(int j=0;j<n;j+=i){
			ll w=1,x,y;
			for(int k=j;k<j+now;k++,w=(w*wn)%mod){
				x=a[k],y=w*a[k+now]%mod;
				a[k]=(x+y)%mod;a[k+now]=((x-y)%mod+mod)%mod;
			}
		}
	}
	if(f==-1){
		ll Inv=fpow(n,mod-2);
		for(int i=0;i<=n;i++) a[i]=(a[i]*Inv)%mod;
		for(int i=1;i<=(n>>1);i++)swap(a[i],a[n-i]);
	}
}

void calc(ll *a,ll *b,int n){
	NTT(a,n,1);NTT(b,n,1);
	for(int i=0;i<=n;i++)a[i]=(a[i]*b[i])%mod;
	NTT(a,n,-1);
}

ll A[M],B[M];
void cdq(int l,int r){
	if(l==r) return;
	int mid=l+r>>1;
	cdq(l,mid);

	int up=r-l-1;
//	for(lg=0,sze=1;sze<=(up<<1);sze<<=1)++lg;--lg;
	for(lg=0,sze=1;sze<=up;sze<<=1)++lg;--lg;//其实只用的到区间长度，实际却有两倍，所以保险用两倍。
	for(int i=0;i<sze;i++) R[i]=(R[i>>1]>>1)|((i&1)<<lg),A[i]=B[i]=0;
//	memset(A,0,sizeof(A));memset(B,0,sizeof(B));
//	fill(A,A+sze,0);fill(B,B+sze,0);//这两种方法赋值fill会比memset快一些（因为fill确定了大小）
	for(int i=l;i<=mid;i++) A[i-l]=f[i];
	for(int i=1;i<=r-l;i++) B[i-1]=g[i];
	calc(A,B,sze);
	
	for(int i=mid+1;i<=r;i++) f[i]=(f[i]+A[i-l-1]%mod)%mod;
	cdq(mid+1,r);
}
int n;
int main(){
	scanf("%d",&n);f[0]=1;
	for(int i=1;i<n;i++)scanf("%lld",&g[i]);
	for(int i=2,up=n<<2;i<=up;i<<=1)inv[i]=fpow(G,(mod-1)/i);
	cdq(0,n-1);
	for(int i=0;i<n;i++) printf("%lld%c",f[i],i==n-1?'\n':' ');
	return 0;
}

---

有没有更加优秀的做法呢？

有，但是只有在取模的意义下，我们只需将式子变形即可。

我们令$F$为$f$的生成函数，$G$为$g$的生成函数。

那么可以得知：

$F(x)=\sum\limits_{i=0}^{\infty}f(i)x^i$

$G(x)=\sum\limits_{i=0}^{\infty}g(i)x^i$

那么求$F$和$G$的卷积就有：

$F(x)\bigotimes G(x)=\sum\limits_{i=0}^{\infty}x^i\sum\limits_{j+k=i}f(j)g(k)$

我们发现$x^i\sum\limits_{j+k=i}f(j)g(k)$是不是和$F(i)$很像，所以我们可以知道那个卷积式可以写成$F(x)-f(0)x^0$，也就是$F(x)-f(0)$，那么我们最后求的答案变为下面这个式子：
$F(x)\bigotimes G(x)\equiv F(x)-f(0)(\mod x^n) \\ F(x)\equiv \frac{f(0)}{1-G(x)}(\mod x^n) \\ F(x)\equiv (1-G(x))^{-1}(\mod x^n)$

那么我们直接使用多项式求逆即可解决，复杂度$O(nlogn)$

---

博主也是才学习，有错误或者疑惑可以提出，大佬带带蒟蒻