学习笔记：ACM算法之FFT（个人理解向）

【音源】【PV確認用】FFT（編集有り）

引言

最近在整理自己学习的算法中发现关于多项式的算法的题已经好久没刷过了。之前也曾学习过此算法，但也仅停留在会用板子上。基于此，我花了数小时的时间重新根据其他大佬的教程自己推导了一遍，然后自己手动实现了一遍FFT的算法。在这里，我也分享一下我自己的一些笔记与推导方法。

FFT概述

FFT(Fast Fourier Transform)，即快速傅里叶变换，是用来加快多项式卷积的计算过程的。想一想两个多项式，如两个3项的多项式：

$$\begin{array}{r}A=1+x+4x^2\\ B=5+x+4x^2\end{array}\text{\hspace{1em}(1)(2)}$$

令 $C=A\times B$，如果我们朴素地计算则为：
$C=1\times (5+x+4x^2)+x\times (5+x+4x^2)+4x^2\times (5+x+4x^2)$

即 $C=5+6x+25x^2+8x^3+16x^4$，需要9次计算（$A$的每一项乘$B$都要三次）。

我们将$A,B$展开为5项(使之与$C$的次数相同)：

$$\begin{array}{r}A=1+x+4x^2+0x^3+0x^4\\ B=5+x+4x^2+0x^3+0x^4\end{array}\text{\hspace{1em}(1)(2)}$$

可推导出$C$的每项系数$c_i$与$A,B$每项系数$a_i,b_i$的关系：
$c_0=a_0\times b_0=5$
$c_1=a_1\times b_0+a_0\times b_1=1+5=6$
$......$
$c_4=a_4\times b_0+a_3\times b_1+a_2\times b_2+a_1\times b_3+a_0\times b_4=0+0+16+0+0=16$

即可推导出：

$$c_i=\sum _{j=0}^i{a}_{i-j}{b}_j\text{\hspace{1em}(3)}$$

刚才我们仅仅算了三项的多项式，如果我们计算$n(n\ge 10^5)$项的多项式，如果还是按照朴素方法计算，时间复杂度为$O(n^2)$，就不能在短时间内计算出来了。这时候，FFT就能派上用场了，它能将原来的时间复杂度$O(n^2)$优化为$O(nlogn)$，也就极大的加快了计算的过程。

一些要用的前置知识

(1)点值表示法

我们考虑线代学习的$(n-1)$次多项式$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$，
众所周知，$n$个互不相同的$x_i$可以唯一确定这个多项式的系数$a_i$。
关于上述的证明，可以考虑将这$n$个等式构成一个非齐次线性方程组，即：

$$\left(\begin{array}{cccccc}1& x_0& x_0^2& \cdots & x_0^{n-2}& x_0^{n-1}\\ 1& x_1& x_1^2& \cdots & x_1^{n-2}& x_1^{n-1}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& x_{n-2}& x_{n-2}^2& \cdots & x_{n-2}^{n-2}& x_{n-2}^{n-1}\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-2}& x_{n-1}^{n-1}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ \cdots \\ a_{n-2}\\ a_{n-1}\end{array}\right)=\left(\begin{array}{c}f(x_0)\\ f(x_1)\\ \cdots \\ f(x_{n-2})\\ f(x_{n-1})\end{array}\right)$$

由于最左端的矩阵为$n\times n$的矩阵，所以可以计算其行列式。又因为其为范德蒙行列式，可求其值为${\prod }_{0\le j<i<n}(x_i-x_j)$。由于$x_i$互不相同，则该行列式的值不为$0$，即该非齐次线性方程组必有唯一解，$n$个互不相同的$x_i$可以唯一确定这个多项式的系数$a_i$。

所以我们可以用这$n$个互不相同的$x_i$组成$n$个点来表示这个多项式，将其写作$((x_0,f(x_0)),(x_1,f(x_1)),\cdots ,(x_{n-1},f(x_{n-1})))$。

(2)复数及单位根

说起复数，相信大家应该都不陌生吧，在我们平常遇到的方程如 $x^2+1=0$，我们容易第一反应这个方程是无解的，这是因为我们第一反应会想到这个方程在实数域上的解。不过，在复数中引入了虚数单位 $i$，且定义 $i^2=-1$，这样就可以得到上述方程其中一个解为$x=i$。

这样，我们就很容易计算两个复数的相乘：

$$\begin{array}{r}A=a+bi\\ B=c+di\end{array}\text{\hspace{1em}(1)(2)}$$

得到：

$$C=A\times B=ac-bd+(bc+ad)i\text{\hspace{1em}(3)}$$

单位根指的是方程 $x^n=1$在复数意义下的解，可以证明，这样的解有$n$个。

我们令 ${\omega }_n=e^{\frac{2\pi i}{n}}$。 这$n$个解为$x={\omega }_n^i(i=0,1,\cdots ,n-1)$，$i$ 表示 ${\omega }_n$ 的幂。从而，我们可以计算得到一下性质：

$$\begin{array}{rl}& x^n={\omega }_n^{ni}=1\\ & {\omega }_n^k={\omega }_{2n}^{2k}\\ & {\omega }_{2n}^{k+n}=-{\omega }_{2n}^k\end{array}\text{\hspace{1em}(1)(2)(3)}$$

在这里简单证明一下：

$(1)$式由欧拉公式$e^{ix}=cosx+isinx$得：
${\omega }_n^{ni}=e^{2\pi i}=cos(2\pi )+isin(2\pi )=1$

$(2)$式有：
${\omega }_{2n}^{2k}=e^{\frac{2\pi i}{2n}\cdot 2k}=e^{\frac{2\pi i}{n}\cdot k}={\omega }_n^k$

$(3)$式由欧拉公式有：
${\omega }_{2n}^{k+n}=e^{\frac{2\pi i}{2n}\cdot (k+n)}=e^{\pi i}\cdot e^{\frac{2\pi i}{2n}\cdot k}=-1\cdot {\omega }_{2n}^k=-{\omega }_{2n}^k$

详细步骤

多项式分组

我们令$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$，根据$x$次数的奇偶将其分为：

$$\begin{array}{r}g(x)=a_0+a_{2}x+a_4{x}^2+\cdots +a_{n-2}{x}^{\frac{n-2}{2}}\\ h(x)=a_1+a_{3}x+a_5{x}^2+\cdots +a_{n-1}{x}^{\frac{n-2}{2}}\end{array}\text{\hspace{1em}(1)(2)}$$

则可以得到：

$$\begin{array}{r}f(x)=g(x^2)+xh(x^2)\end{array}\text{\hspace{1em}(3)}$$

对于$(3)$式，我们把${\omega }_n^k和{\omega }_n^{k+\frac{n}{2}}$代入，计算得：

$$\begin{array}{rl}f({\omega }_n^k)& =g(({\omega }_n^k{)}^2)+{\omega }_n^{k}h(({\omega }_n^k{)}^2)\\ & =g({\omega }_n^{2k})+{\omega }_n^{k}h({\omega }_n^{2k})\\ & =g({\omega }_{\frac{n}{2}}^k)+{\omega }_n^{k}h({\omega }_{\frac{n}{2}}^k)\end{array}$$

$$\begin{array}{r}f({\omega }_n^{k+\frac{n}{2}})=f(-{\omega }_n^k)=g({\omega }_{\frac{n}{2}}^k)-{\omega }_n^{k}h({\omega }_{\frac{n}{2}}^k)\end{array}\text{\hspace{1em}(1)}$$

所以在我们知道$g({\omega }_{\frac{n}{2}}^k)$和$h({\omega }_{\frac{n}{2}}^k)$后，就可以同时计算出$f({\omega }_n^k)$和$f({\omega }_n^{k+\frac{n}{2}})$的值了。

FFT正变换

对于FFT正变换，得到的目标为$f({\omega }_n^i)(i=0,1,\cdots ,n-1)$这$n$个值，如果我们按照$f(x)$多项式本身，把 ${\omega }_n^i$ 一个一个代入计算，时间复杂度还是$O(n^2)$。这个时候，我们就要利用多项式分组中得到的性质，自底向上递推出这$n$个值了。

在这里举个例子，令 $n=8,f(x)=a_0+a_{1}x+\cdots +a_7{x}^7$，对于$f(x)$，我们先将其分组：

$(1)$第一次分组后，
$$\begin{array}{r}g(x)=a_0+a_{2}x+a_4{x}^2+a_6{x}^3\\ h(x)=a_1+a_{3}x+a_5{x}^2+a_7{x}^3\end{array}\text{\hspace{1em}(1)(2)}$$

$(2)$我们将上述$g(x),h(x)$，继续分组，令$g(x)$得到的分组为$g_1^{\prime }(x),h_1^{\prime }(x)$，$h(x)$得到的分组为$g_2^{\prime }(x),h_2^{\prime }(x)$，即可得到：
$$\begin{array}{r}{g}_1^{\prime }(x)=a_0+a_{4}x\\ h_1^{\prime }(x)=a_1+a_{5}x\\ g_2^{\prime }(x)=a_2+a_{6}x\\ h_2^{\prime }(x)=a_3+a_{7}x\end{array}\text{\hspace{1em}(1)(2)(3)(4)}$$

$(3)$以此类推直到每个函数都只有一项：
$$\begin{array}{r}{g}_1^{\prime \prime }(x)=a_0\\ h_1^{\prime \prime }(x)=a_4\\ g_2^{\prime \prime }(x)=a_1\\ h_2^{\prime \prime }(x)=a_5\\ g_3^{\prime \prime }(x)=a_2\\ h_3^{\prime \prime }(x)=a_6\\ g_4^{\prime \prime }(x)=a_3\\ h_4^{\prime \prime }(x)=a_7\end{array}\text{\hspace{1em}(1)(2)(3)(4)(5)(6)(7)(8)}$$

分组后就可以变换了：

这里先把前面所述的两个计算公式再重复一遍，以便阅读：
$f({\omega }_n^k)=g({\omega }_{\frac{n}{2}}^k)+{\omega }_n^{k}h({\omega }_{\frac{n}{2}}^k)$
$f({\omega }_n^{k+\frac{n}{2}})=g({\omega }_{\frac{n}{2}}^k)-{\omega }_n^{k}h({\omega }_{\frac{n}{2}}^k)$

$(1)$刚开始先从${\omega }_1^0$开始计算，按照下标相同的为一组代入公式可得到：
$f^{\prime }({\omega }_2^0)=g^{\prime \prime }({\omega }_1^0)+{\omega }_2^0{h}^{\prime \prime }({\omega }_1^0)$
$f^{\prime }({\omega }_2^1)=g^{\prime \prime }({\omega }_1^0)-{\omega }_2^0{h}^{\prime \prime }({\omega }_1^0)$

总共会计算$8$次，得到的$f^{\prime }$作为下一次计算所使用的$g^{\prime },h^{\prime }$，即：
$g_1^{\prime }({\omega }_2^0),h_1^{\prime }({\omega }_2^0),g_2^{\prime }({\omega }_2^0),h_2^{\prime }({\omega }_2^0),g_1^{\prime }({\omega }_2^1),h_1^{\prime }({\omega }_2^1),g_2^{\prime }({\omega }_2^1),h_2^{\prime }({\omega }_2^1)$

$(2)$将得到的$g^{\prime },h^{\prime }$按照$(1)$中方法再次计算(即计算${\omega }_2^0,{\omega }_2^1$)，得到：
$g({\omega }_4^0),g({\omega }_4^1),g({\omega }_4^2),g({\omega }_4^3),h({\omega }_4^0),h({\omega }_4^1),h({\omega }_4^2),h({\omega }_4^3)$
总共也是计算$8$次。

$(3)$最后再进行一次，即可得到$f({\omega }_n^i)(i=0,1,\cdots ,n-1)$这$n$个值。

我们可以从上述推导过程得出，对于$n$项多项式，我们需要将其规整成$2^m$的形式，即 $n\le 2^m<2n$，对于$n\ne 2^m$的$n$项的多项式来说，则需在其高次项补零，即：
$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}+0x^n+0x^{n+1}+\cdots +0x^{2^m-1}$

同时，对于上述自底向上的计算，可以证明每次递推的计算次数为$2^m$次(或$n$次)，共进行$logn$次递推。所以对于FFT的正变换，时间复杂度为$O(nlogn)$。

FFT逆变换

我们回顾一下点值表示法对应的非齐次线性方程组，将系数$a_i$作为未知数，代入$f({\omega }_n^i)(i=0,1,\cdots ,n-1)$，可以得到：

$$\left(\begin{array}{cccccc}1& 1& 1& \cdots & 1& 1\\ 1& {\omega }_n^1& {\omega }_n^2& \cdots & {\omega }_n^{n-2}& {\omega }_n^{n-1}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& {\omega }_n^{n-2}& {\omega }_n^{2(n-2)}& \cdots & {\omega }_n^{(n-2{)}^2}& {\omega }_n^{(n-2)(n-1)}\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-2)(n-1)}& {\omega }_n^{(n-1{)}^2}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ \cdots \\ a_{n-2}\\ a_{n-1}\end{array}\right)=\left(\begin{array}{c}f(1)\\ f({\omega }_n^1)\\ \cdots \\ f({\omega }_n^{n-2})\\ f({\omega }_n^{n-1})\end{array}\right)$$

由范德蒙矩阵的逆，可得到系数$a_i$为：

$$\left(\begin{array}{c}{a}_0\\ a_1\\ \cdots \\ a_{n-2}\\ a_{n-1}\end{array}\right)=\frac{1}{n}\left(\begin{array}{cccccc}1& 1& 1& \cdots & 1& 1\\ 1& {\omega }_n^{-1}& {\omega }_n^{-2}& \cdots & {\omega }_n^{-(n-2)}& {\omega }_n^{-(n-1)}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& {\omega }_n^{-(n-2)}& {\omega }_n^{-2(n-2)}& \cdots & {\omega }_n^{-(n-2{)}^2}& {\omega }_n^{-(n-2)(n-1)}\\ 1& {\omega }_n^{-(n-1)}& {\omega }_n^{-2(n-1)}& \cdots & {\omega }_n^{-(n-2)(n-1)}& {\omega }_n^{-(n-1{)}^2}\end{array}\right)\left(\begin{array}{c}f(1)\\ f({\omega }_n^1)\\ \cdots \\ f({\omega }_n^{n-2})\\ f({\omega }_n^{n-1})\end{array}\right)$$

即可得到$a_i=\frac{1}{n}(f(1)+f({\omega }_n^1){\omega }_n^{-i}+f({\omega }_n^2)({\omega }_n^{-i}{)}^2+\cdots +f({\omega }_n^{n-1})({\omega }_n^{-i}{)}^{n-1})$，可以看成$a_i=\frac{1}{n}f({\omega }_n^{-i})$的形式，同正变换一样，再做一次FFT即可。

个人用板子

本人使用的是递归版本的FFT，来源：快速傅里叶变换

#include <math.h>
#include <iostream>
#include <complex>

using namespace std;

typedef complex<double> Comp;

const Comp I(0, 1); //Comp I(a, b) = a + b * i
const int maxn = 1 << 20; //根据实际题目大小自行调整
const double M_PI_ = acos(-1.0);

Comp temp[maxn], f[maxn];

void DFT(Comp* f, int n, int rev) { //0 ~ n - 1 共n项 , rev = 1, DFT; rev = -1, IDFT, 计算f(wnk)
	if (n == 1) { //常数项
		return;
	}
	for (int i = 0; i < n; i++) {
		temp[i] = f[i];
	}
	for (int i = 0; i < n; i++) {
		if (i & 1) {
			f[i / 2 + n / 2] = temp[i]; // 奇数位置范围：n / 2 ~ n - 1;
		}
		else {
			f[i / 2] = temp[i]; // 偶数位置范围：0 ~ n / 2 - 1;
		}
	}
	Comp* g = f, * h = f + n / 2;
	Comp cur(1, 0);
	Comp step = exp(I * (2 * M_PI_ / n * rev));
	DFT(g, n / 2, rev);
	DFT(h, n / 2, rev);
	for (int k = 0; k < n / 2; ++k) {
		temp[k] = g[k] + cur * h[k];
		temp[k + n / 2] = g[k] - cur * h[k];
		cur *= step;
	}
	for (int i = 0; i < n; ++i) {
		f[i] = temp[i];
	}
}

两道例题

P3803 【模板】多项式乘法（FFT）

思路：分别对多项式$F(x)$和$G(x)$进行FFT正变换，得到的值分别相乘，作为逆变换的$f$，最后逆变换即可。

#include <math.h>
#include <iostream>
#include <complex>

using namespace std;

typedef complex<double> Comp;

const Comp I(0, 1); //Comp I(a, b) = a + b * i
const int maxn = 1 << 24;
const double M_PI_ = 3.1415926535;

Comp temp[maxn], f[maxn], g[maxn];

void DFT(Comp* f, int n, int rev) { //0 ~ n - 1 共n项 , rev = 1, DFT; rev = -1, IDFT, 计算f(wnk)
	if (n == 1) { //常数项
		return;
	}
	for (int i = 0; i < n; i++) {
		temp[i] = f[i];
	}
	for (int i = 0; i < n; i++) {
		if (i & 1) {
			f[i / 2 + n / 2] = temp[i]; // 奇数位置范围：n / 2 ~ n - 1;
		}
		else {
			f[i / 2] = temp[i]; // 偶数位置范围：0 ~ n / 2 - 1;
		}
	}
	Comp* g = f, * h = f + n / 2;
	Comp cur(1, 0);
	Comp step = exp(I * (2 * M_PI_ / n * rev));
	DFT(g, n / 2, rev);
	DFT(h, n / 2, rev);
	for (int k = 0; k < n / 2; ++k) {
		temp[k] = g[k] + cur * h[k];
		temp[k + n / 2] = g[k] - cur * h[k];
		cur *= step;
	}
	for (int i = 0; i < n; ++i) {
		f[i] = temp[i];
	}
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(0);
	int n, m;
	cin >> n >> m;
	int size = 1;
	int max_ = n + m;
	while (size <= max_) {
		size <<= 1;
	}
	for (int i = 0; i <= n; i++) {
		int a;
		cin >> a;
		f[i].real(a);
	}
	for (int i = 0; i <= m; i++) {
		int a;
		cin >> a;
		g[i].real(a);
	}
	DFT(f, size, 1);
	DFT(g, size, 1);
	for (int i = 0; i < size; i++) {
		double freal = f[i].real(), greal = g[i].real(), fimag = f[i].imag(), gimag = g[i].imag();
		f[i].real(freal * greal - fimag * gimag);
		f[i].imag(freal * gimag + fimag * greal);
	}
	DFT(f, size, -1);
	for (int i = 0; i <= n + m; i++) {
		cout << (int)(f[i].real() / size + 0.5) << " ";
	}
	return 0;
}

P3338 [ZJOI2014]力

思路：把公式转换为卷积的形式，最终可得到$E_j={\sum }_{i=0}^{j}f(i)g(j-i)-{\sum }_{i=0}^{n-j+1}g(i)h(n-j+1-i)$，其中$h$为$f$的$n$到$0$的序列(即反向序列)。

#include <math.h>
#include <iostream>
#include <complex>

using namespace std;

typedef complex<double> Comp;

const Comp I(0, 1); //Comp I(a, b) = a + b * i
const int maxn = 1 << 20; //根据实际题目大小自行调整
const double M_PI_ = acos(-1.0);

Comp temp[maxn], f[maxn], g[maxn], h[maxn];

void DFT(Comp* f, int n, int rev) { //0 ~ n - 1 共n项 , rev = 1, DFT; rev = -1, IDFT, 计算f(wnk)
	if (n == 1) { //常数项
		return;
	}
	for (int i = 0; i < n; i++) {
		temp[i] = f[i];
	}
	for (int i = 0; i < n; i++) {
		if (i & 1) {
			f[i / 2 + n / 2] = temp[i]; // 奇数位置范围：n / 2 ~ n - 1;
		}
		else {
			f[i / 2] = temp[i]; // 偶数位置范围：0 ~ n / 2 - 1;
		}
	}
	Comp* g = f, * h = f + n / 2;
	Comp cur(1, 0);
	Comp step = exp(I * (2 * M_PI_ / n * rev));
	DFT(g, n / 2, rev);
	DFT(h, n / 2, rev);
	for (int k = 0; k < n / 2; ++k) {
		temp[k] = g[k] + cur * h[k];
		temp[k + n / 2] = g[k] - cur * h[k];
		cur *= step;
	}
	for (int i = 0; i < n; ++i) {
		f[i] = temp[i];
	}
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(0);
	int n;
	cin >> n;
	int size = 1;
	int max_ = 2 * n;
	while (size <= max_) {
		size <<= 1;
	}
	for (int i = 1; i <= n; i++) {
		double a;
		cin >> a;
		f[i].real(a);
		h[n - i + 1].real(a);
		g[i].real(1.0 / (double)i / (double)i);
	}
	DFT(f, size, 1);
	DFT(g, size, 1);
	DFT(h, size, 1);
	for (int i = 0; i <= size; i++) {
		double freal = f[i].real(), greal = g[i].real(), hreal = h[i].real();
		double fimag = f[i].imag(), gimag = g[i].imag(), himag = h[i].imag();
		f[i].real(freal * greal - fimag * gimag);
		f[i].imag(freal * gimag + fimag * greal);
		g[i].real(hreal * greal - himag * gimag);
		g[i].imag(hreal * gimag + himag * greal);
	}
	DFT(f, size, -1);
	DFT(g, size, -1);
	for (int i = 1; i <= n; i++) {
		printf("%.6lf\n", (f[i].real() - g[n - i + 1].real()) / size);
	}
	return 0;
}

总结

这应该是我第一次正式记录自己学习的算法，目的还是为了加深对这个变换思想的理解。此博客仅为自己的学习记录，可能有一些理解不同或者错误的地方，也请各位大佬多多指教啦~
顺便浇浇FFT的紫谱和白谱喵