快速傅里叶变换(FFT)的分治实现

本文作者为 JustinRochester。

目录地址

上一篇

下一篇


快速傅里叶变换(FFT)的分治实现

引入

根据 第一篇 的定义,当我们拥有两个度数分别不超过 \(n,m\) (两者同阶)的多项式 \(\displaystyle F(x)=\sum_{i=0}^n f_ix^i, G(x)=\sum_{i=0}^m g_ix^i\) 时,我们可以计算两者的多项式卷积:

\(\displaystyle H(x)=(F\cdot G)(x)=\sum_{i=0}^{n+m} x^i\sum_{p+q=i}f_pg_q\)

在这种情况下,我们计算两个多项式的卷积,即是通过长为 \(n+1, m+1\) 的数组 \(f, g\) 求解长为 \(n+m+1\) 的数组 \(h\) ,使得 \(\displaystyle h_k=\sum_{i+j=k}f_ig_j\)

在朴素的算法下,我们求解的复杂度是 \(O(n^2)\) 的。

但根据 第二篇 的知识,如果能求出 \(H(x)\) 上的 \(n+m+1\) 个点,我们通过拉格朗日插值可以唯一地求解一个不超过 \(n+m\) 次的多项式 \(H(x)\)

虽然暂时这个方法的复杂度是一致的 \(O(n^2)\) 且常数更大,但它提供了一个更好的思路:我们能不能把多项式用一定数量的点来进行表示?


多项式的点值表示

我们规定,我们用 \(n+1\) 个点(即 \(x\) 点的值与对应的 \(y\) 值)来表示一个度数不超过 \(n\) 的多项式 \(\displaystyle F(x)=\sum_{i=0}^n f_ix^i\) 。这称为 \(F(x)\) 的点值表示法。

朴素上,我们可以在 \(O(n^2)\) 的时间内,通过秦九韶算法求解这些点的值。

当我们需要计算 \(H(x)\) 的乘积时,我们需要用 \(n+m+1\) 个点表示这个不超过 \(n+m\) 次的多项式;而恰好,由于 \(F(x), G(x)\) 的次数都不超过 \(n+m\) ,因此也可以用 \(n+m+1\) 个点进行表示。

而当我们求解 \(H(x)\) 的点值表示结果时,假设我们已经通过上述算法求解出了 \(F(x), G(x)\) 的点值表示 \((x_1, F(x_1)), (x_1, G(x_1)), (x_2, F(x_2)), (x_2, G(x_2)), \cdots\) ,我们可以直接在 \(O(n)\) 的时间内求解 \(H(x)\) 的点值表示法:

\((x_i, H(x_i))=(x_i, F(x_i)\cdot G(x_i))\) ,只要对应采样点的 \(y\) 值相乘即可。

于是我们得到了一个不太一样的多项式加减乘运算的方法:

\(F(x), G(x)\stackrel{\text{秦九韶算法采样}}{\longrightarrow}F(x), G(x)\text{ 的点值表示}\stackrel{\text{点值乘积}}{\longrightarrow}H(x)\text{ 的点值表示}\stackrel{\text{拉格朗日插值}}{\longrightarrow}H(x)\)

这个算法的第一步和第三步是 \(O(n^2)\) 的,而第二步是 \(O(n)\) 的。


快速傅里叶变换

很显然,这个算法的第二步并没有什么优化的空间了。于是,我们考虑如何优化第一步和第三步(下一篇的内容)。

对于第一步的采样,如果我们采样的点之间没有相互的关联,我们只能通过秦九韶算法硬采样。那如果我们代入一些有关联的点,比如单位根呢?

我们假设现在对一个不超过 \(n-1\) 次的多项式 \(\displaystyle F(x)=\sum_{i=0}^{n-1} f_ix^i\) 代入 \(n\)\(n\) 次单位根的各幂次。则 \(F(x)\) 的第 \(i\) 个点值表示法的结果为:\((\omega_n^i, F(\omega_n^i))\)

为了方便,我们不妨设 \(n\) 为偶数(不足的可以前补 \(0\))。我们令 \(\displaystyle F_L(x)=\sum_{i=0}^{{n\over 2}-1}f_{2i}x^i, F_R(x)=\sum_{i=0}^{{n\over 2}-1}f_{2i+1}x^i\)

这样一来,我们可以很容易的发现,\(F(x)=F_L(x^2)+xF_R(x^2)\) ,实际上 \(F_L(x),F_R(x)\) 就是将数组 \(F(x)\) 进行奇偶分列。

我们再次代入采样点,会发现 \(\displaystyle F(\omega_n^i)=F_L(\omega_n^{2i})+\omega_n^iF_R(\omega_n^{2i})=F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i)\)

我们不妨假设 \(i<{n\over 2}\) ,则再检查 \(\omega_n^{i+{n\over 2}}\) 的结果,就会同样有:

\(\displaystyle F(\omega_n^{i+{n\over 2}})=F_L(\omega_{n\over 2}^{i+{n\over 2}})+\omega_n^{i+{n\over 2}}F_R(\omega_{n\over 2}^{i+{n\over 2}})\)

这里同样根据 第四篇 的内容,我们知道:\(\displaystyle \omega_{n\over 2}^{i+{n\over 2}}=\omega_{n\over 2}^i, \omega_n^{i+{n\over 2}}=-\omega_n^i\)

于是我们代入化简就会得到:\(\displaystyle F(\omega_n^{i+{n\over 2}})=F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i)\)

对比两式,我们归纳得到:当 \(i<{n\over 2}\) 时:

\(\begin{cases} \begin{aligned} F(\omega_n^i)&=&F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i)\\ F(\omega_n^{i+{n\over 2}})&=&F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i)\\ \end{aligned} \end{cases}\)

因此我们发现,求解小于 \(n\) 次的多项式采样点 \(F(\omega_n^i)\) 的问题,就简化为了如何求解两个小于 \({n\over 2}\) 次的多项式采样点 \(F_L(\omega_{n\over 2}^i), F_R(\omega_{n\over 2}^i)\) 的问题。

而后者一旦求解完毕,我们可以在 \(O(n)\) 的时间内,合并得到所有的 \(F(\omega_n^i)\)

而对于后者的求解,我们会发现,它的求解结构和原问题一模一样!我们只需要继续递归就可以了。

递归的边界在于多项式减小到不超过 \(1\) 次时,我们已经有 \(F_i(x)=f_i\) 了。此时的递归可以直接返回结果,\(F_i(\omega_1^1)=f_i\)

于是时间复杂度为 \(\displaystyle T(n)=2T({n\over 2})+\Theta(n)=\Theta(n\log n)\) 。这个复杂度的证明我们放在后续。

这种实现方法一般叫做 DIT-FFT ,即基于时间抽取的 FFT。另一个常见的实现方法我们会在后续介绍,是 DIF-FFT,即基于频率抽取的 FFT。


快速傅里叶变换的分治实现

当然,在递归求解的过程中,也可能出现次数变为奇数的问题。这个问题可以在递归的过程中,像上文那样补充一项的 \(0\) 使得其变成偶数。但我们一般的实现过程中,为了避免这种麻烦,我们一般都是直接将 \(n\) 补充为 \(2^k\)

这个方法的另一个好处是:我们可以预处理所有需要用到的单位根。例如 \(\displaystyle \omega_{n\over 8}^i=\omega_{2^{k-3}}^i=\omega_{2^k}^{8i}=\omega_n^{8i}\) ,我们只需要修改步长即可。

若采用每次升为最近偶数的方法,当多项式次数为 \(17\) 时,我们补充到 \(18\) 并获得了 \(\omega_{18}^0,\cdots, \omega_{18}^{17}\)

而第一次递归分解为了两个次数为 \(9\) 的多项式,我们又需要补充到 \(10\) ,然后需要重新计算 \(\omega_{10}^0,\cdots, \omega_{10}^9\) ,因为 \(\omega_{10}^i\) 不能全部在 \(\omega_{18}^j\) 中找到;

继续第二次递归后,得到两个次数不超过 \(5\) 的,又需要补充到 \(6\) ,然后重新计算 \(\omega_6^0,\cdots, \omega_6^5\) ......

依次类推,每次都需要重新计算单位根,加大了运算量。

因此我们将快速傅里叶变换的分治实现总结如下:

输入:
一个度数小于 \(n\) 的多项式 \(\displaystyle F(x)=\sum_{i=0}^{n-1}f_ix^i\)

预处理:
\(n\) 补充为 \(2^k\) ,预处理 \(\omega_n^0, \cdots, \omega_n^{n-1}\)

过程:

  1. \(n=1\) 则直接返回;
  2. \(F(x)\) 奇偶分列为 \(F_L(x), F_R(x)\)
  3. 递归求解 \(F_L(x), F_R(x)\) 的结果;
  4. \(i<{n\over 2}\) ,计算 \(\begin{cases} \begin{aligned} F(\omega_n^i)&=&F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i)\\ F(\omega_n^{i+{n\over 2}})&=&F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i)\\ \end{aligned} \end{cases}\)
  5. 返回计算结果;

快速傅里叶变换的分治实现-参考代码

int N;
vir w[MAXN];//单位根
inline void FFT(vir *f, int n) {
	static vir tmp[MAXN];
	if(n==1) return ;
	
	//奇偶分列
	for(int i=0; i<n; ++i) tmp[i]=f[i];
	vir *fl=f, *fr=f+n/2;
	for(int i=0, j=0; i<n; i+=2, ++j) fl[j]=tmp[i];
	for(int i=1, j=0; i<n; i+=2, ++j) fr[j]=tmp[i];
	
	//递归求解
	FFT(fl, n/2); FFT(fr, n/2);
	
	//O(n) 合并
	vir *o=w;
	int pace = N/n;
	for(int i=0; i<n/2; ++i) {
		tmp[i] = fl[i] + *o * fr[i]; 
		tmp[i + n/2] = fl[i] - *o * fr[i];
		o += pace;
	}
	for(int i=0; i<n; ++i) f[i]=tmp[i];
}
inline void mul(vir *f, vir *g, int n, int m) {
	for(N=1; N<n+m-1; N<<=1);//求最小的 2 的幂次
	vir omega = get_omega(N);
	
	w[0]=1;
	for(int i=1; i<N; ++i) w[i]=w[i-1]*omega;
	FFT(f, N); FFT(g, N);
	for(int i=0; i<N; ++i) f[i]*=g[i];
	IFFT(f, N);
}

附录:FFT 复杂度的简单求解

根据 FFT 的复杂度计算式:\(\displaystyle T(n)=2T({n\over 2})+\Theta(n)\)

我们不妨假设,当 \(n\leq N\) 时,\(T(n)\leq C_1\) ,以及 \(n\geq N\)\(\displaystyle T(n)\leq 2T({n\over 2})+C_2n\)

\(n=N\cdot 2^k\) 时,有:

\(\begin{aligned} &T(n)\\ =&T(N\cdot 2^k)\\ \leq&2T(N\cdot 2^{k-1})+C_2\cdot N\cdot 2^k\\ \leq&4T(N\cdot 2^{k-2})+2\cdot C_2\cdot N\cdot 2^{k-1}+C_2\cdot N\cdot 2^k\\ \leq&4T(N\cdot 2^{k-2})+2\cdot C_2\cdot N\cdot 2^k\\ \leq&\cdots\\ \leq&2^kT(N)+k\cdot C_2\cdot N\cdot 2^k\\ \leq&2^k\cdot C_1+k\cdot C_2\cdot N\cdot 2^k\\ =&{C_1\over N}\cdot n+C_2(n\cdot \log_2{n\over N})\\ =&C_2\cdot n\log_2 n+({C_1\over N}-C_2\cdot \log_2 N)\cdot n\\ \end{aligned}\)

同理,不妨假设 \(n\geq N_1\)\(\displaystyle T(n)\geq 2T({n\over 2})+C_3n\) 可得到 \(\displaystyle T(n)\geq C_3\cdot n\log_2 n - C_3\cdot \log_2 N\cdot n\)

因此,\(\exists n_0, c_d=C_3-1, c_u=C_2+1\text{ s.t. }n\geq n_0\to c_d\cdot n\log_2 n\leq T(n)\leq c_u\cdot n\log_2 n\)\(T(n)=\Theta(n\log n)\)

否则,我们令 \(n_0\leq N\cdot 2^{k-1}<n<N\cdot 2^k\) ,有:

\(n<N\cdot 2^k\)\(T(n)\leq T(N\cdot 2^k)\leq c_u\cdot k\cdot 2^k< c_u\cdot {n\over N}\cdot \log_2 {n\over N}={c_u\over N}\cdot n\log_2 n - {c_u\over N}\log_2 N\cdot n\leq {c_u\over N}\cdot n\log_2 n\leq c_u\cdot n\log_2 n\)

另一方向 \(n>N\cdot 2^{k-1}\) 同理可得 \(T(n)\geq c_d\cdot n\log_2 n\)

于是,此时也有 \(T(n)=\Theta(n\log n)\)

综上,该算法的复杂度是 \(\Theta(n\log n)\) 的。

posted @ 2023-02-06 10:30  JustinRochester  阅读(342)  评论(0编辑  收藏  举报