学习总结-快速傅里叶变换（FastFourierTransform）

\(update\ on\ 2020.9.11:\)​增加了“实战演练”（放题解的地方）

\(update\ on\ 2021.7.31:\)增加了分治法

（〇）前言

快速傅里叶变换（Fast Fourier Transform, FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法。FFT会通过把DFT矩阵分解为稀疏（大多为零）因子之积来快速计算此类变换。因此，它能够将计算DFT的复杂度从只用DFT定义计算需要的复杂度\(\Theta(N^2)\)，降低到\(\Theta(N\log N)\)，其中\(N\)为数据大小。

快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及，但早在1805年就已推导出来。1994年美国数学家吉尔伯特·斯特朗把FFT描述为“我们一生中最重要的数值算法”，它还被IEEE科学与工程计算期刊列入20世纪十大算法。

不要管上面那段话。作为OIer，我们的关注点应该只有两个：

1. 离散傅里叶变换（DFT）是什么

2. 时间复杂度\(\Theta(N^2)\to\Theta(N\log N)\)

现在就由小编带领大家快速入门FFT。

前置知识：实数的运算（逃

（一）从多项式到卷积

1.1多项式的系数表示法与点值表示法

多项式（Polynomial）是代数学中的基础概念，是由称为未知数的变量和称为系数的常数通过有限次加减法、乘法以及自然数幂次的乘方运算得到的代数表达式。

多项式是整式的一种。未知数只有一个的多项式称为一元多项式。下文讨论的多项式指一元多项式。

给出一个多项式：

\[A(x)=\sum\limits_{i=0}^{N-1}a_ix^i \]

不难发现，当\(a_0,a_1,...,a_{N-1}\)确定了以后，多项式\(A(x)\)也就确定了。因此，我们可以用\((a_0,a_1,...,a_{N-1})\)来表示一个多项式。这就是多项式的系数表示法。

将\(N\)个互不相同的\(x\)，分别为\((x_0,x_1,...,x_{N-1})\)（称为插值节点），带入到多项式\(A(x)\)中，就可以得到\(N\)组值，分别为\((A(x_0),A(x_1),...,A(x_{N-1}))\)，或者表示为\((y_0,y_1,...,y_{N-1})\)。这样一来，多项式也随之确定。这就是多项式的点值表示法。

理解点值表示法可以参考\(N\)元一次方程组。把\(N\)个\(x\)不同的取值带进去，就得到了\(N\)个\(N\)元一次方程。分别为：

\(\begin{bmatrix}x_0^0a_0,&x_0^1a_1,&...,&x_0^{N-2}a_{N-2},&x_0^{N-1}a_{N-1}\\x_1^0a_0&x_1^1a_1,&...,&x_1^{N-2}a_{N-2},&x_1^{N-1}a_{N-1}\\\vdots&\vdots&\vdots&\vdots&\vdots\\x_{N-1}^0a_{0}&x_{N-1}^1a_1&...,&x_{N-1}^{N-2}a_{N-2}&x_{N-1}^{N-1}a_{N-1}\end{bmatrix}=\begin{bmatrix}y_0\\y_1\\\vdots\\y_{N-1}\end{bmatrix}\)

于是\(a_0,a_1,...,a_{N-1}\)都唯一确定了，即多项式也唯一确定了。解出上述方程就能求得系数表示法下的多项式了。

在这里有一个重要的结论：系数表示法和点值表示法是可以相互转化的。

1.2多项式乘法与卷积

现在有两个多项式\(A(x)=\sum\limits_{i=0}^{N-1}a_ix^i,B(x)=\sum\limits_{i=0}^{M-1}b_ix^i\)

如果我们把上述两个多项式乘起来，那么我们就会得到一个\(N+M-2\)次多项式。设这个多项式为\(C(x)\)

那么，

\[C(x)=A(x)\times B(x)=\sum\limits_{k=0}^{N+M-2}\sum\limits_{i+j=k}a_ib_jx^k \]

解释一下：由于\(C(x)\)是一个\(N+M-2\)次多项式，那么他的所有项的次数就为\(0,1,2...N+M-4,N+M-3,N+M-2\)。因此我们只需要枚举每一项的次数，算出该项的系数就行了。

举个例子:

\(c_5x^5=a_0x^0b_5x^5+a_1x^1b_4x^4+...+a_5x^5b_0x^0=(a_0b_5+a_1b_4+...+a_5b_0)x^5\)

以上就是多项式乘法。实际上就是将两个多项式的项两两相乘，简记为\(C=A*B\)。

接下来就是乘法的本质——卷积。

设\(c_k\)为\(C\)的第\(k\)项的系数。\(a_k,b_k\)同理。

由多项式乘法可以得到：

\[c_k=\sum\limits_{i+j=k}a_ib_j \]

形如\(c_k=\sum\limits_{i\oplus j=k}a_ib_j\)的式子称为卷积，其中\(\oplus\)是某种运算。多项式乘法实际上就是加法卷积。

卷积还可以写作函数的形式，就像\(C=A*B\)。

介绍另一种卷积：狄利克雷卷积（在我的学习总结-莫比乌斯反演中有描述）

\[(f*g)(n)=\sum\limits_{d|N}f(d)g(\dfrac{N}{d}) \]

狄利克雷卷积是乘法卷积（\(d|N\)实际上就是\(i\times j=N,i=d,j=\dfrac{N}{d}\)）

1.3求解多项式乘法

求多项式\(C=A*B\)。

方法很简单，直接按照定义计算就行了。给出代码：

for(int i=0; i<N; i++)
	for(int j=0; j<M; j++)
		c[i+j] += a[i]*b[j];

就是上面这个入门难度的东西啦。时间复杂度：\(\Theta(N^2)\)

然而，讲这么简单的东西干什么。是因为——这个东东可以用FFT优化到\(\Theta(N\log N)\)！没错，一开始我听说这个东西的时候也很震惊：这玩意儿还能优化？

欲后事如何，且听下回分解。

（二）复数的引入

写本章时学校突然断电，然后gg

2.1奇妙的分治

珂学家们经过研究发现，当用点值表示多项式的时候，有一些好用的性质：

对于一组插值节点\((x_0, x_1, x_2,...,x_{N-1}),\)

多项式\(A\)的值为\((A(x_0),A(x_1),A(x_2),...,A(x_{N-1}))\)

多项式\(B\)的值为\((B(x_0),B(x_1),B(x_2),...,B(x_{N-1}))\)

多项式\(C=A*B\)的值为\((C(x_0),C(x_1),C(x_2),...,C(x_{N-1}))\)

那么，有\(C(x_0)=A(x_0)\times B(x_0),C(x_1)=A(x_1)\times B(x_1),...,C(x_{N-1})=A(x_{N-1})\times B(x_{N-1}))\)

也就是说，在点值表示法下，我们可以\(\Theta(N)\)求出两个多项式的乘积。

然而，多项式一般是用系数表示法表示的。问题就转化为了：如何快速将多项式从系数表示法转化为点值表示法，在求得乘积后，如何快速将点值表示法转化为系数表示法。

我们知道，系数表示法转点值表示法要通过插入\(N\)个不同的插值节点。而每次插值的时间复杂度是\(\Theta(N)\)，总时间复杂度就是\(\Theta(N^2)\)，我们貌似又回到了起点。

此时，珂学家又提出了一个算法：分治。即将一个多项式\(A\)按照项的奇偶性分为两个子问题求解。于是就有了（假设多项式的项数\(N\)为\(2^k\)）：

\(A(x)=(a_0x^0+a_2x^2+a_4x^4+...+a_{N-2}x^{N-2})+(a_1x^1+a_3x^3+a_5x^5+...+a_{N-1}x^{N-1})\)

\(=(a_0x^0+a_2x^2+a_4x^4+...+a_{N-2}x^{N-2})+x(a_1x^0+a_3x^2+a_5x^4+...+a_{N-1}x^{N-2})\)

设：

\[A_1=\sum\limits_{i=0}^{\frac{N-2}{2}}a_{2i}x^i,A_2(x)=\sum\limits_{i=0}^{\frac{N-2}{2}}a_{2i+1}x^i \]

可以得到：

\[A(x)=A_1(x^2)+xA_2(x^2) \]

这时候，

\[C(x)=(A_1(x^2) + xA_2(x^2))(B_1(x^2)+xB_2(x^2)) \]

拆开看一看：

\[C(x)=A_1(x^2)B_1(x^2)+x(A_1(x^2)B_2(x^2)+A_2(x^2)B_1(x^2))+x^2A_2(x^2)B_2(x^2) \]

如果直接这样计算的话，需要计算4次多项式乘法。

时间复杂度为 \(T(N)=\Theta(N)+4T(\frac{N}{2})\)​，大概是 \(\Theta(N^2)\)​​。

然而还可以继续化简：

\[C(x)=(1-x)A_1(x^2)B_1(x^2)+x(A_1(x^2)+A_2(x^2))(B_1(x^2)+B_2(x^2))+(x^2-x)A_2(x^2)B_2(x^2) \]

乍一看貌似更复杂了，事实上不难发现，计算这个式子只需要三次多项式乘法。

时间复杂度为 \(T(N)=\Theta(N)+3T(\frac{N}{2})\)，大约是 \(\Theta(N^{1.47})\)，十分优秀。

2.2为什么要引入复数

当我们在算\(x\)​​的时候，用到了\(x^2\)​​的插值结果。如果\(x^2\)​​也是插值节点，是否可以直接调用？

我们可以想象一下我们需要怎样的插值节点：

• 如果\(x\)存在，那么\(x^2\)也存在。
• 插值节点的个数要等于\(N\)

可以从复平面内找到这\(N\)个特殊的插值节点。

以上是引入复数的原因。

2.3复数的基本认识

（1）复数

复数，为实数的延伸，它使任一多项式方程都有根。

复数当中有个“虚数单位”\(i\)，它是\(-1\)的一个平方根，即\(i^2=-1\)。任一复数都可表达为\(a+bi\)，其中\(a\)及\(b\)皆为实数，分别称为复数之“实部”和“虚部”。

每个复数都对应了平面上的一个向量\((a,b)\)，这个平面叫做复平面，它是由水平的实轴与垂直的虚轴建立起来的复数的几何表示。

故每一个复数唯一对应了一个复平面上的向量，每一个复平面上的向量也唯一对应了一个复数。其中\(0\)既被认为是实数，也被认为是虚数。

复数\(z\)的模长\(|z|\)定义为\(z\)在复平面到原点的距离，\(|z|=\sqrt{a^2+b^2}\)。幅角\(\theta\)为实轴的正半轴正方向（逆时针）旋转到\(z\)的有向角度。

由于虚数无法比较大小。复数之间的大小关系只存在等于与不等于两种关系，两个复数相等当且仅当实部虚部对应相等。对于虚部为\(0\)的复数之间是可以比较大小的，相当于实数之间的比较。

复数之间的运算满足结合律，交换律和分配律。由此定义复数之间的运算法则：

\[(a+bi)+(c+di)=(a+c)+(b+d)i \]

\[(a+bi)-(c+di)=(a-c)+(b-d)i \]

\[(a+bi)\cdot(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i \]

\[\dfrac{a+bi}{c+di}=\dfrac{(a+bi)\cdot(c-di)}{(c+di)\cdot(c-di)}=\dfrac{ac+bd}{c^2+d^2}+\dfrac{(bc-ad)i}{c^2+d^2} \]

复数运算的加法满足平行四边形法则，乘法满足幅角相加，模长相乘。

对于一个复数\(z=a+bi\)，它的共轭复数是\(z'=a-bi\)，\(z'\)称为\(z\)的复共轭。

共轭复数有一些性质:

\(z\cdot z'=a^2+b^2\)

\(|z|=|z'|\)

复数类的代码实现（仅供参考）：

struct Complex{
	double r, i;
	Complex() {}
	Complex(double r, double i): r(r), i(i) {}
	inline Complex operator +(const Complex& x) const {
		return Complex(r + x.r, i + x.i);
	}
	inline Complex operator -(const Complex& x) const {
		return Complex(r - x.r, i - x.i);
	}
	inline Complex operator *(const Complex& x) const {
		return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline void operator *=(const Complex& x) {
		*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline Complex conj(){
		return Complex(r, -i);
	}
};

（2）单位根

下图是复数中的单位圆。

根据欧拉幅角公式，\(e^{i\varphi}=\cos\varphi+i\sin\varphi\)

我们发现单位圆上的复数（设幅角为\(\varphi\)）可以表示为\(e^{i\varphi}\)。

将单位圆等分成\(N\)个部分（以单位圆与实轴正半轴的交点一个等分点），以原点为起点，圆的这\(N\)个\(N\)等分点为终点，作出\(N\)个向量。这\(N\)个向量就是复数的单位根啦。

其中幅角为正且最小的向量称之为\(N\)次单位向量，记为\(\omega_n^1\)。

以此类推，\(N\)个向量分别为\(\omega_n^1,\omega_n^2,\omega_n^3,...,\omega_n^n\)，它们可以用复数之间的乘法得来\(\omega_n^k=\omega_n^{k-1}\cdot\omega_n^1(2\leq k\leq n)\)

显然\(\omega_n^0=\omega_n^n=1\)

\(\omega_n^k\)也可以写作\(e^{i\cdot \frac{2k\pi}{n}}\)

所以\(\omega_n^k=\cos(\dfrac{2k\pi}{n})+i\sin(\dfrac{2k\pi}{n})\)

\(\omega_n^k\)有两个重要的性质：

1. 性质1：折半引理

\[\omega_{2n}^{2k}=\omega_{n}^{k} \]

由几何意义可以知道，这两者表示的向量终点是一样的。

同样可以用三角恒等变换证明。

2. 性质2：消去引理

\[\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k} \]

由几何意义可以知道，两者表示的向量模长相等，方向相反。

也可以用三角恒等变换证明。

这里给出求\(\omega_n^k\)的代码实现：

omega[i] = Complex(cos(2.0*PI/N*i), sin(2.0*PI/N*i));

（三）快速傅里叶变换（Fast Fourier Transform）

3.1离散傅里叶变换（Discrete Fourier Transform）

（1）理论分析

步入正题。如何使用复数加速插值。

相信经过上一张的讨论，大家心里都已经清楚2.1中提到的这\(N\)个特殊的复数是什么了。这\(N\)个复数正是复数的单位根。

不难发现，复数的\(N\)个单位根恰好满足了我们“理想”中插值节点的所有条件。温习一下：

• 如果\(x\)存在，那么\(x^2\)也存在。
• 插值节点的个数要等于\(N\)

接下来我们要将这\(N\)个单位根带入2.1中推导出的公式。公式如下：

\[A(x)=A_1(x^2)+xA_2(x^2) \]

分为两种情况：

设\(0\leq k<\dfrac{n}{2}\)，则大于等于\(\dfrac{n}{2}\)的数可以表示为\(k+\dfrac{n}{2}\)

1. 对于\(0\leq k<\dfrac{n}{2}\)的部分

\(A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^{k}A_2(\omega_n^{2k})\)

由折半引理得：\(\omega_n^{2k}=\omega_\frac{n}{2}^k\)

\(\therefore A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_\frac{n}{2}^k)\)

2. 对于\(\dfrac{n}{2}\leq k+\dfrac{n}{2}<n\)的部分

\(A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})\)

由消去引理得：\(\omega_n^{2k+n}=\omega_n^{2k},\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)

由折半引理得：\(\omega_n^{2k}=\omega_\frac{n}{2}^{k}\)

\(\therefore A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_\frac{n}{2}^{k})-\omega_n^kA_2(\omega_\frac{n}{2}^k)\)

把上述两个式子合起来：

\[\begin{cases}A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_\frac{n}{2}^k)\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_\frac{n}{2}^{k})-\omega_n^kA_2(\omega_\frac{n}{2}^k)\end{cases}(0\leq k\leq\dfrac{n}{2}) \]

依照上面那条式子，我们就能以\(\Theta(N)\)的时间复杂度带入每一个根\(\omega_n^k\)，然后转化为子问题\(A_1,A_2\)求解。\(A_1,A_2\)的规模都是原问题的一半，因此，时间复杂度为：

\(T(n)=2T(\dfrac{N}{2})+\Theta(N)=\Theta(N\log_2 N)\)

以上是离散傅里叶变换的理论，可以快速将系数表示法转化为点值表示法。

（2）程序实现

上述操作可以通过递归实现。下面给出C++代码：

代码中的部分细节：

\(\operatorname{Complex}\)是自己封装的复数类（用\(\operatorname{STL}\)也可以，不过比较慢）

\(\operatorname{omega}\)是预处理出来的\(\omega\)值，\(\operatorname{omega}[k]\)存的是\(\omega_n^k\)的值。

由于在递归的过程中\(n\)会改变，设\(n\)改变后变为\(len\)，那么根据折半引理，\(\omega_{len}^k=\omega_{n}^{\frac{n}{len}k}\)

代码应该还是比较清晰的。

void transform(Complex *a, const int& len, const Complex* omega, const int& N){
	if(len == 1) return ;
	int mid = len>>1;
	for(int i=0; i<len; i++) save[i] = a[i];
	Complex *al = a, *ar = a + mid;
	for(int i=0; i<len/2; i++) al[i] = save[i<<1], ar[i] = save[i<<1|1];
	transform(al, len>>1, omega, N), transform(ar, len>>1, omega, N);
	for(int i=0; i<len/2; i++){
		save[i] = al[i] + ar[i]*omega[N/len*i];
		save[i+len/2] = al[i] - ar[i]*omega[N/len*i];
	}
	for(int i=0; i<len; i++) a[i] = save[i];
}

以上。

3.2离散傅里叶反变换（Inverse Discrete Fourier Transform）

（1）理论分析

有了DFT（离散傅里叶变换）还不够，我们还需要将点值表示法转化为系数表示法。这时候，我们就需要用到IDFT（离散傅里叶反变换）。

也就是我们要从\((A(x_0),A(x_1),A(x_2),...,A(x_{n-1}))\)推导出系数\((a_0,a_1,a_2,...,a_{n-1})\)

进入莫比乌斯反演模式

先将\((a_0,a_1,a_2,...,a_{n-1})\)做DFT，得到\((b_0,b_1,b_2,...,b_{n-1})\)

写成多项式形式：\(B(x)=\sum\limits_{i=0}^{n-1}b_ix^i\)，有：

\[b_k=\sum\limits_{i=0}^{n-1}a_i\cdot(\omega_n^k)^i \]

将\(x=\omega_n^{-k}\)带入\(B(x)\)，得到一组点值\((c_0,c_1,c_2, ...,c_{n-1})\)

就得到了这样的式子：

\[c_k=\sum\limits_{i=0}^{n-1}b_i\cdot(\omega_i^{-k})^i \]

然后将两式合并：

\(c_k=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j\cdot(\omega_n^i)^j\cdot(\omega_i^{-k})^i\)

\(=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{-k})^i\cdot(\omega_n^{i})^j\)

\(=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}\omega_n^{ij-ki}\)

\(=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{i})^{j-k}\)

设\(\delta =j-k,S(j,k)=\sum\limits_{i=0}^{n-1}(\omega_n^i)^\delta=\sum\limits_{i=0}^{n-1}(\omega_n^\delta)^i\)

1. 当\(\omega_n^\delta=1\)时，\(\delta=0,j=k,S(j,k)=n\)

2. 当\(\omega_n^\delta\not=1\)时，根据等比数列求和公式，\(S(j,k)=\dfrac{\omega_n^0[(\omega_n^\delta)^{n}-1]}{\omega_n^\delta-1}=\dfrac{1\times[(\omega_n^n)^{\delta}-1]}{\omega_n^\delta-1}=\dfrac{1-1}{\omega_n^\delta-1}=0\)

\(\therefore c_k=\sum\limits_{j=0}^{n-1}a_j\cdot n\cdot[j=k]=a_k\cdot n\)

神奇的事情发生了...我们得到了一个十分简洁美丽的公式：

\[a_k=\dfrac{c_k}{n} \]

其中\(c_k\)是插值节点\(\omega_n^{-k}\)代入到\(B(x)\)后得到的值。

以上就是IDFT的理论。

简单地说，DFT就是用\(\omega_n^k\)插值系数表达得到点值表达，IDFT就是用\(\omega_n^{-k}\)插值系数表达得到点值表达。

（2）代码实现

与DFT的实现除了插值变为了\(\omega_n^{-k}\)以外完全相同（当然得到的结果还要除以\(n\)）。

给出调用的代码：

\(\operatorname{omegaInverse}\)数组指的是\(\omega_n^{-k}\)，直接用共轭复数的求法就可以求出。（详细见2.2（1））

void idft(Complex *a, const int& N){
    transform(a, N, omegaInverse, N);
    for(int i=0; i<N; i++) a[i]/=N;
}

以上。

3.3FFT加速多项式乘法

所谓FFT，其实就是将DFT和IDFT合起来。下面是加速多项式乘法的流程：

1. 使用DFT，将\(A(x),B(x)\)转为点值表示。复杂度\(\Theta(N\log N)\)
2. 在点值表示下将两个多项式相乘，得到\(C(x)\)的点值表示。复杂度\(\Theta(N)\)
3. 使用IDFT，将\(C(x)\)转为系数表示。复杂度\(\Theta(N\log N)\)

FFT总之间复杂度\(\Theta(N\log N)\)，而且常数较大。

下面是递归FFT的代码实现（上文大部分已经解释过了）：

网上大部分题解的递归FFT都无法通过模板多项式乘法

因为已经预处理出了\(\omega_n^k\)和\(\omega_n^{-k}\)，下面的程序常数较小，实测可以通过。

#include<bits/stdc++.h>
using namespace std;
const double PI = 3.1415926535;
const int maxn = 4e6 + 5;
inline double read(){
	int w = 0, f = 1; char ch = getchar();
	while(ch < '0' or ch > '9') {if(ch == '-') f = -f; ch = getchar();}
	while(ch >= '0' and ch <= '9') w = w*10 + ch - '0', ch = getchar();
	return (double)w*f;
}
struct Complex{
	double r, i;
	Complex() {}
	Complex(double r, double i): r(r), i(i) {}
	inline Complex operator +(const Complex& x) const {
		return Complex(r + x.r, i + x.i);
	}
	inline Complex operator -(const Complex& x) const {
		return Complex(r - x.r, i - x.i);
	}
	inline Complex operator *(const Complex& x) const {
		return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline void operator *=(const Complex& x) {
		*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline Complex conj(){
		return Complex(r, -i);
	}
}A[maxn], B[maxn], C[maxn];
struct FastFourierTransform{
	Complex omega[maxn], omegaInverse[maxn], save[maxn];
	void init(const int& N){
		for(int i=0; i<N; i++){
			omega[i] = Complex(cos(2.0*PI/N*i), sin(2.0*PI/N*i));
			omegaInverse[i] = omega[i].conj();
		}
	}
	void transform(Complex *a, const int& len, const Complex* omega, const int& N){
		if(len == 1) return ;
		int mid = len>>1;
		for(int i=0; i<len; i++) save[i] = a[i];
		Complex *al = a, *ar = a + mid;
		for(int i=0; i<len/2; i++) al[i] = save[i<<1], ar[i] = save[i<<1|1];
		transform(al, len>>1, omega, N), transform(ar, len>>1, omega, N);
		for(int i=0; i<len/2; i++){
			save[i] = al[i] + ar[i]*omega[N/len*i];
			save[i+len/2] = al[i] - ar[i]*omega[N/len*i];
		}
		for(int i=0; i<len; i++) a[i] = save[i];
	}
	void dft(Complex *a, const int& N){
		transform(a, N, omega, N);
	}
	void idft(Complex *a, const int& N){
		transform(a, N, omegaInverse, N);
	}
}fft;
int main(){
	int N, M;
	N = read(), M = read();
	for(int i=0; i<=N; i++) A[i].r = read();
	for(int j=0; j<=M; j++) B[j].r = read();
	int limit = 1;
	while(limit <= N+M) limit <<= 1;
	fft.init(limit);
	fft.dft(A, limit); fft.dft(B, limit);
	for(int i=0; i<limit; i++) C[i] = A[i]*B[i];
	fft.idft(C, limit);
	for(int i=0; i<=N+M; i++) printf("%d ", (int)(C[i].r/limit + 0.5)); 
	return 0;
}

以上。其实到这里已经可以了但我们要追求极致的效率。

3.4精细FFT

为什么叫精细FFT呢因为它很精细呀。为了追求极致的效率，我们选择对上文的递归实现FFT进行卡常优化。

（1）蝴蝶操作

操作的名称由来是因为这个操作很像蝴蝶飞舞（可以想象这个操作是多么的秀）。

观察下图，显示了递归的过程。将序号转为二进制。

不难发现，递归最底层的数所在的位置恰好是值的二进制翻转。

比如说，\(6\)的二进制是\(110\)，而它在递归底层的位置是\(3\)，二进制是\(011\)，恰好是\(6\)的二进制的翻转。

据此，我们可以将最底层的所有值的位置确定，然后逐层向上合并，这样能够大大减小常数。这就是蝴蝶操作。

问题来了，如何求某个数的二进制翻转。

用类似于数位DP的方式，设\(\operatorname{reserve(x)}\)表示\(x\)的二进制翻转。

1. 当\(x\)的二进制最后一位为\(1\)时，

reserve[x]=reserve[x>>1]>>1|(N>>1);

2. 当\(x\)的二进制最后一位为\(0\)时，

reserve[x]=reserve[x>>1]>>1;

合起来就是：

reserve[x] = (reserve[x>>1]>>1)|((x&1)?N>>1:0);

解释一下，如果我们要翻转一个\(x\)位的二进制数，就要把最后一位拿出来，把剩下的\(x-1\)位翻转，再把拿出来的最后一位放到第一位。如下：

			@****1
拿走1，变成		@****
翻转，变成		****@
把1放回，变成	        1****@

代码上面已经有了。

（2）三步变两步

观察复数的乘法：

\((a+bi)^2=a^2+2abi-b^2\)

设多项式\(P(x)=A(x)+iB(x)\)

那么\(P(x)^2=(A(x)^2+iB(x)^2)=A(x)^2-B(x)^2+2A(x)B(x)i\)

我们发现，多项式\(C=A*B\)的两倍恰好就是\(P(x)^2\)的虚部！

因此我们只需要把\(A(x)\)当做\(P(x)\)的实部，\(B(x)\)当做\(P(x)\)的虚部，然后求\(P(x)^2\)，就能将\(C(x)\)求出啦。

这样就只做了一次DFT，一次IDFT，由三步变为了两步。

奉上本人最精细的FFT:

#include<bits/stdc++.h>
using namespace std;
const int maxn = 4e6 + 5;
inline double read(){
	int w = 0, f = 1; char ch = getchar();
	while(ch < '0' or ch > '9') {if(ch == '-') f = -f; ch = getchar();}
	while(ch >= '0' and ch <= '9') w = w*10 + ch - '0', ch = getchar();
	return (double)w*f;
}
double PI = 3.1415926535;
struct Complex{
	double r, i;
	Complex() {}
	Complex(double r, double i) : r(r), i(i) {}
	inline void real(const double& x) {r = x;}
	inline double real() {return r;}
	inline Complex operator + (const Complex& rhs) const {
		return Complex(r + rhs.r, i + rhs.i);
	}
	inline Complex operator - (const Complex& rhs) const {
		return Complex(r - rhs.r, i - rhs.i);
	}
	inline Complex operator * (const Complex& rhs) const {
		return Complex(r*rhs.r - i*rhs.i, r*rhs.i + i*rhs.r);
	}
	inline void operator /= (const double &x) {
		r /= x, i /= x;
	}
	inline void operator *= (const Complex& rhs) {
		*this = Complex(r*rhs.r - i*rhs.i, r*rhs.i + i*rhs.r);
	}
	inline void operator += (const Complex& rhs) {
		r += rhs.r, i += rhs.i;
	}
	inline Complex conj() {
		return Complex(r, -i);
	}
};

struct FastFourierTransform{
	Complex omega[maxn], omegaInverse[maxn]; int reserve[maxn];
	void init(const int& N){
		for(int i=0; i<N; i++){
			omega[i] = Complex(cos(2*PI/N*i), sin(2*PI/N*i));
			omegaInverse[i] = omega[i].conj();
			reserve[i] = (reserve[i>>1]>>1)|((i&1)?N>>1:0);
		}
	}
	void transform(Complex *a, const int& N, const Complex* omega){
		for(int i=0; i<N; i++) if(reserve[i] > i) swap(a[reserve[i]], a[i]);
		for(int l=2; l<=N; l<<=1){
			int m = l/2;
			for(int j = 0; j < N; j += l){
				for(int i=0; i<m; i++) {
					Complex t = omega[N/l*i] * a[j+m+i];
					a[j+m+i] = a[j + i] - t;
					a[j + i] = a[j + i] + t;
				}
			}
		}
	}
	void dft(Complex *a, const int& N){
		transform(a, N, omega);
	}
	void idft(Complex *a, const int& N){
		transform(a, N, omegaInverse);
	}
}fft;
Complex F[maxn];
int main(){
	int N, M;
	N = read(), M = read();
	for(int i=0; i<=N; i++) F[i].r = read();
	for(int j=0; j<=M; j++) F[j].i = read();
	int limit = 1;
	while(limit <= N+M) limit <<= 1;
	fft.init(limit);
	fft.dft(F, limit);
	for(int i=0; i<limit; i++) F[i] = F[i]*F[i];
	fft.idft(F, limit);
	for(int i=0; i<=N+M; i++) printf("%d ", (int)(F[i].i/limit/2 + 0.5));
	return 0;
}

以上。全剧终。

（四）实战演练

经过了社会的毒打我发现了解FFT和使用FFT是两回事，于是便有了这一章...

1.1[ZJOI2014]力

题目是求\(E_i\)的值：

\(F_j=\sum\limits_{i=1}^{j-1}\dfrac{q_i\times q_j}{(i-j)^2}-\sum\limits_{i=j+1}^{n}\dfrac{q_i\times q_j}{(i-j)^2},E_i=\dfrac{F_i}{q_i}\)

第一眼看上去：这跟FFT有什么关系？——哦，原来要化简。

那就化简吧：

\(E_j=\dfrac{F_j}{q_j}=\sum\limits_{i=1}^{j-1}\dfrac{q_i}{(j-i)^2}-\sum\limits_{i=1}^{n-j}\dfrac{q_{i+j}}{i^2}\)

观察前面的式子\(\sum\limits_{i=1}^{j-1}\dfrac{q_i}{(j-i)^2}\)，如果设\(f(i)=q_i,g(i)=\dfrac{1}{i^2}\)，那么就变成了：

\(\sum\limits_{i=1}^{j-1}f(i)g(j-i)\)，写好看一点：

\[\sum\limits_{i=0}^{j}f(i)g(j-i) \]

为了不影响结果，我们要令\(g(0)=0,f(j)=0\)。

上面的式子其实已经变成了卷积的形式。FFT的作用是求多项式乘法，多项式乘法事实上就是卷积。

也就是说，上面那个东西可以用FFT算出来（与多项式乘法无异）。

但是后面还有一个式子：\(\sum\limits_{i=1}^{n-j}\dfrac{q_{i+j}}{i^2}\)
，让它顺眼一点：

\[\sum\limits_{i=0}^{n-j}f(i+j)g(i) \]

同样为了保证不影响结果，我们要令\(f(j)=0,g(0)=0\)。

这个式子就有点神奇了，因为它没有\(i+j=k\)的形式（不是卷积的形式）。

然而这是有套路的。虽然\(i+j+i\not=n-j\)，但是我们发现\(i+j-i=j\)。所以我们要想办法让\(i+j+i\to i+j-i\)。

两种方法：

1. 令\(g'(n-i)=g(i)\)，原式\(=\sum\limits_{i=0}^{n-j}f(i+j)g'(n-i)=\sum\limits_{i=0}^{n+j}f(i+j)g'(n-i)\)

2. 令\(f'(n-i-j)=f(i+j)\)，原式\(=\sum\limits_{i=0}^{n-j}f'(n-i-j)g(i)\)

其实就是将一个数组翻转，将差变为和。以上两种方法任君选择。

方法一需要注意的是，为了不影响结果，要保证\(\sum\limits_{i=n-j+1}^{n+j}f(i+j)g'(n-i)=0\)

最后，分别用FFT算出和式的结果后，作差即可。

两种方法在细节处理上略有差别。

附上代码（两种方法的代码都有）：

（方法一）

#include<bits/stdc++.h>
using namespace std;
const int maxn = 4e5 + 5;
const double PI = acos(-1);
struct Complex{
	double r, i;
	Complex () {}
	Complex (double r, double i): r(r), i(i) {}
	inline Complex operator +(const Complex& x) const{
		return Complex(r + x.r, i + x.i);
	}
	inline Complex operator -(const Complex& x) const{
		return Complex(r - x.r, i - x.i);
	}
	inline Complex operator *(const Complex& x) const{
		return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline void operator +=(const Complex& x){
		*this = Complex(r + x.r, i + x.i);
	}
	inline void operator *=(const Complex& x){
		*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline Complex conj(){
		return Complex(r, -i);
	}
};
struct FastFourierTransform{
	Complex omega[maxn], omegaInverse[maxn]; int reserve[maxn];
	void init(const int& N){
		for(int i=0; i<N; i++){
			omega[i] = Complex(cos(2*PI/N*i), sin(2*PI/N*i));
			omegaInverse[i] = omega[i].conj();
			reserve[i] = (reserve[i>>1]>>1)|((i&1)?N>>1:0);
		}
	}
	void transform(Complex *a, const int& N, const Complex* omega){
		for(int i=0; i<N; i++) if(i > reserve[i]) swap(a[i], a[reserve[i]]);
		for(int l=2; l<=N; l<<=1){
			int m = l>>1;
			for(Complex *p=a; p != a+N; p += l)
				for(int i=0; i<m; i++){
					Complex t = p[i+m]*omega[N/l*i];
					p[i+m] = p[i] - t;
					p[i] = p[i] + t;
				}
		}
	}
	void dft(Complex *a, const int& N){
		transform(a, N, omega);
	}
	void idft(Complex *a, const int& N){
		transform(a, N, omegaInverse);
	}
}fft;
Complex q[maxn], g1[maxn], g2[maxn];
int main(){
	int N;
	scanf("%d", &N);
	for(int i=1; i<=N; i++) scanf("%lf", &q[i].r);
	for(int i=1; i<=N; i++) g1[i].r = 1.0/(double)i/i, g2[N-i].r = g1[i].r;
	int limit = 1;
	while(limit <= N<<1) limit <<= 1;
	fft.init(limit);
	fft.dft(g1, limit);
	fft.dft(g2, limit);
	fft.dft(q, limit);
	for(int i=0; i<limit; i++) g1[i] = q[i]*g1[i];
	for(int i=0; i<limit; i++) g2[i] = q[i]*g2[i];
	fft.idft(g1, limit);
	fft.idft(g2, limit);
	for(int i=1; i<=N; i++){
		printf("%.3lf\n", (g1[i].r - g2[N+i].r)/(double)limit);
	}
	return 0;
}

（方法二）

#include<bits/stdc++.h>
using namespace std;
const int maxn = 4e5 + 5;
const double PI = acos(-1);
struct Complex{
	double r, i;
	Complex () {}
	Complex (double r, double i): r(r), i(i) {}
	inline Complex operator +(const Complex& x) const{
		return Complex(r + x.r, i + x.i);
	}
	inline Complex operator -(const Complex& x) const{
		return Complex(r - x.r, i - x.i);
	}
	inline Complex operator *(const Complex& x) const{
		return Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline void operator +=(const Complex& x){
		*this = Complex(r + x.r, i + x.i);
	}
	inline void operator *=(const Complex& x){
		*this = Complex(r*x.r - i*x.i, r*x.i + i*x.r);
	}
	inline Complex conj(){
		return Complex(r, -i);
	}
};
struct FastFourierTransform{
	Complex omega[maxn], omegaInverse[maxn]; int reserve[maxn];
	void init(const int& N){
		for(int i=0; i<N; i++){
			omega[i] = Complex(cos(2*PI/N*i), sin(2*PI/N*i));
			omegaInverse[i] = omega[i].conj();
			reserve[i] = (reserve[i>>1]>>1)|((i&1)?N>>1:0);
		}
	}
	void transform(Complex *a, const int& N, const Complex* omega){
		for(int i=0; i<N; i++) if(i > reserve[i]) swap(a[i], a[reserve[i]]);
		for(int l=2; l<=N; l<<=1){
			int m = l>>1;
			for(Complex *p=a; p != a+N; p += l)
				for(int i=0; i<m; i++){
					Complex t = p[i+m]*omega[N/l*i];
					p[i+m] = p[i] - t;
					p[i] = p[i] + t;
				}
		}
	}
	void dft(Complex *a, const int& N){
		transform(a, N, omega);
	}
	void idft(Complex *a, const int& N){
		transform(a, N, omegaInverse);
	}
}fft;
Complex q1[maxn], q2[maxn], g[maxn];
int main(){
	int N;
	scanf("%d", &N);
	for(int i=1; i<=N; i++) scanf("%lf", &q1[i].r), q2[N-i].r = q1[i].r;
	for(int i=1; i<=N; i++) g[i].r = 1.0/((double)i*i);
	int limit = 1;
	while(limit <= N<<1) limit <<= 1;
	fft.init(limit);
	fft.dft(g, limit);
	fft.dft(q1, limit);
	fft.dft(q2, limit);
	for(int i=0; i<limit; i++) q1[i] = q1[i]*g[i];
	for(int i=0; i<limit; i++) q2[i] = q2[i]*g[i];
	fft.idft(q1, limit);
	fft.idft(q2, limit);
	for(int i=1; i<=N; i++){
		printf("%.3lf\n", (q1[i].r - q2[N-i].r)/(double)limit);
	}
	return 0;
}

未完待续...

（五）参考文献

1. 知乎 作者：白空谷 文章：一小时学会快速傅里叶变换（Fast Fourier Transform）

2. command_block 的博客 文章：傅里叶变换(FFT)学习笔记

3. wikipedia维基百科-快速傅里叶变换

4. wikipedia维基百科-复数

5. wikipedia维基百科-欧拉公式

上述1.和2.都是很好的文章，可以加深对FFT的理解。