简单理解 FFT

现在有两个序列 \(f_{0\cdots n-1}\) 与 \(g_{0\cdots m-1}\)，我们需要计算 \(h_{0\cdots n+m-2}\) 满足：

\[h_i=\sum_{j=0}^if_jg_{i-j} \]

要求在 \(O((n+m)\log (n+m))\) 的复杂度内完成计算。FFT 就是干这个事情的。

一般来说，我们会先计算 \(f,g\) 的 DFT，即

\[F_i=\sum_{j=0}^{n-1}f_j\omega_n^{ij} \]

其中 \(\omega_{n}^k=\cos \frac{2k\pi}{n}+\text{i}\sin \frac{2k\pi }{n}\) 为 \(n\) 次单位根。它是方程 \(x^n=1\) 的 \(n\) 个复数根。

比较直观的理解方式是把 \(\omega_n^{0\cdots n-1}\) 理解为在单位圆上等距离地转一圈。

那么我们就认为 \(F\) 是 \(f\) 的 DFT。由 \(f\) 求 \(F\) 叫 DFT，由 \(F\) 求 \(f\) 叫 IDFT。显然两个是一一对应的。

• 事实上如果你把 \(f\) 看作多项式 \(f(x)=\sum f_ix^i\) 那么 \(f\) 的 DFT 就是它在 \(w^{0\cdots n-1}\) 处的点值。

FFT 干的事情就是，它分别给 \(f,g\) 取了这样一组点值，然后分别相乘，那么由于 \(h(x)=f(x)\times g(x)\)，相乘后得到的点值就是 \(h\) 的点值，只需做 IDFT 即可得到 \(h\)。

显然关键在于怎么求这组点值。（以及怎么通过点值快速还原）

首先我们一般会给多项式补 \(0\)，找到最小的 \(2^k\ge n+m\)，把多项式都补到 \(2^k\) 位，方便分治。

然后我们进行分治。这个分治比较特殊，对于要做 DFT 的序列 \(f_{0\cdots 2^k-1}\)，它会按照奇偶性分类，也就是把 \(f_0,f_2,f_4,\cdots f_{2^{k}-2}\) 分到一组，\(f_1,f_3,\cdots f_{2^k-1}\) 分到一组。把它们看做两个长度为 \(2^{k-1}\) 的序列计算 DFT，我们考虑怎么通过这两个东西算出 \(f_{0\cdots 2^k-1}\) 的 DFT。

• 递归边界是只有一个数时，它的 DFT 是它本身。这相当于常数多项式。

然后我们推一个奥妙重重的式子：设 \(A(x)=\sum f_ix^i,A_0(x)=\sum f_{2i}x^i,A_1(x)=\sum f_{2i+1}x^i\)

那么有 \(A(x)=A_0(x^2)+xA_1(x^2)\)。如果我们想算 \(A(\omega_{2^k}^i)\)，那么

\[A(\omega_{2^k}^i)=A_0(\omega_{2^k}^{2i})+\omega_{2^k}^iA_1(\omega_{2^k}^{2i})=A_0(\omega_{2^{k-1}}^i)+\omega_{2^k}^iA_1(\omega_{2^{k-1}}^i)\\ A(\omega_{2^k}^{i+2^{k-1}})=A_0(\omega_{2^k}^{2(i+2^{k-1})})+\omega_{2^k}^{i+2^{k-1}}A_1(\omega_{2^k}^{2({i+2^{k-1}})})=A_0(\omega_{2^{k-1}}^i)-\omega_{2^k}^iA_1(\omega_{2^{k-1}}^i) \]

上面两式中均有 \(i<2^{k-1}\)。

其中用到了两条单位根的性质：\(\omega_{n}^k=\omega_{nd}^{kd},\omega_{2m}^k+\omega_{2m}^{k+m}=0\)。读者不难自证。

乍一看很神秘，貌似通过一个玄之又玄的手段计算出了 \(f_{0\cdots 2^{k}-1}\)。但你仔细想一下发现，单位根就是在单位圆上等距地跳 \(n\) 步，以 \(F_1=\sum_{i=0}^{n-1}\omega_n^if_i\) 为例，你发现它其实是把 \(f_i\) 与跳第 \(i\) 步到达的点值对上，然后两两相乘，再全部加起来。

那我们按照奇偶性分类之后，跳的距离恰好是原来的两倍，我们只会跳到 \(\omega_{n}^{0,2,\cdots,n}\)。而刚好分过去的 \(f\) 也是隔一位分一个，于是偶数那边恰好有 \(f_0\) 对 \(\omega_n^0\)，\(f_2\) 对 \(\omega_n^2\)，等等。奇数那边几乎同理，不过需要同时再跳一步，因此第 \(i\) 项需要在奇数那边乘上 \(\omega_{n}^i\)。

那么 IDFT 怎么算呢，我告诉你有这个式子：

\[b_k=\sum_{i=0}^{n-1}a_i\omega_n^{ki}\iff a_k=\dfrac{1}{n}\sum_{i=0}^{n-1}b_i\omega_n^{-ki} \]

你发现这两个式子的区别无非是把 \(\omega_n\) 换成了 \(\omega_n^{-1}\)，因此只需要稍微修改一下实现就行了

然后实现上直接递归空间复杂度 \(O(n\log n)\) 常数又太大，考虑自底向上计算。这个分治树比较特别，我们就考虑把分治树重构成区间的形式，方便自底向上计算。发现只需要把每个位置的二进制表示逆序一下即可。

计算的时候对于区间 \([l,m],[m+1,r]\)，根据推的式子发现只需要按照 \(i=0,1,\cdots,m-l\) 的顺序，每次同时计算 \(F_{l+i},F_{m+1+i}\) 即可。这样就做到了线性空间。

然后你发现这个利用的性质无非是 \(\omega_n\) 的那两条，那么我告诉你 \(\bmod n\) 意义下如果存在原根 \(g\)，那么 \(g^{0,1,\cdots n-1}\) 同样恰好构成循环。因此在模意义下我们也可以做 FFT，这个叫 NTT，通常模数是 \(998244353\)，原根 \(g=3\)，且满足 \(2^{23}\mid 998244353-1\)。

void NTT(int *A,int tag){
	for(int i=0;i<N;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int i=1;i<=L;i++){
		int m=1<<i,wn=ksm((tag==1)?g:ig,(mod-1)/m);
		omega[0]=1;for(int l=1;l<(m>>1);l++)omega[l]=omega[l-1]*wn%mod;
		for(int j=0;j<N;j+=m){
			for(int l=0;l<(m>>1);l++){
				int x=A[j+l],y=A[j+l+(m>>1)]*omega[l]%mod;
				A[j+l]=(x+y)%mod,A[j+l+(m>>1)]=(x-y+mod)%mod;
			}
		}
	}
	for(int i=0;i<N;i++)A[i]=(A[i]%mod+mod)%mod;
	if(tag==-1){int iv=inv(N);for(int i=0;i<N;i++)A[i]=A[i]*iv%mod;}
}