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

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

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

（〇）前言

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

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

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

离散傅里叶变换（DFT）是什么
时间复杂度(Theta(N^2) oTheta(Nlog N))

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

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

（一）从多项式到卷积

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

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

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

给出一个多项式：

[A(x)=sumlimits_{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)元一次方程。分别为：

(egin{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}=egin{bmatrix}y_0\y_1\vdots\y_{N-1}end{bmatrix})

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

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

1.2多项式乘法与卷积

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

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

那么，

[C(x)=A(x) imes B(x)=sumlimits_{k=0}^{N+M-2}sumlimits_{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=sumlimits_{i+j=k}a_ib_j ]

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

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

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

[(f*g)(n)=sumlimits_{d|N}f(d)g(dfrac{N}{d}) ]

狄利克雷卷积是乘法卷积（(d|N)实际上就是(i imes 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(Nlog 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) imes B(x_0),C(x_1)=A(x_1) imes B(x_1),...,C(x_{N-1})=A(x_{N-1}) imes 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=sumlimits_{i=0}^{frac{N-2}{2}}a_{2i}x^i,A_2(x)=sumlimits_{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})。幅角( heta)为实轴的正半轴正方向（逆时针）旋转到(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)的复共轭。

共轭复数有一些性质:

(zcdot 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^{ivarphi}=cosvarphi+isinvarphi)

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

将单位圆等分成(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}cdotomega_n^1(2leq kleq n))

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

(omega_n^k)也可以写作(e^{icdot frac{2kpi}{n}})

所以(omega_n^k=cos(dfrac{2kpi}{n})+isin(dfrac{2kpi}{n}))

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

性质1：折半引理

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

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

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

性质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) ]

分为两种情况：

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

对于(0leq 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)

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

对于(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})

( herefore A(omega_n^{k+frac{n}{2}})=A_1(omega_frac{n}{2}^{k})-omega_n^kA_2(omega_frac{n}{2}^k))

把上述两个式子合起来：

[egin{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}(0leq kleqdfrac{n}{2}) ]

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

(T(n)=2T(dfrac{N}{2})+Theta(N)=Theta(Nlog_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)=sumlimits_{i=0}^{n-1}b_ix^i)，有：

[b_k=sumlimits_{i=0}^{n-1}a_icdot(omega_n^k)^i ]

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

就得到了这样的式子：

[c_k=sumlimits_{i=0}^{n-1}b_icdot(omega_i^{-k})^i ]

然后将两式合并：

(c_k=sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}a_jcdot(omega_n^i)^jcdot(omega_i^{-k})^i)

(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}(omega_n^{-k})^icdot(omega_n^{i})^j)

(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}omega_n^{ij-ki})

(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}(omega_n^{i})^{j-k})

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

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

( herefore c_k=sumlimits_{j=0}^{n-1}a_jcdot ncdot[j=k]=a_kcdot 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合起来。下面是加速多项式乘法的流程：

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

FFT总之间复杂度(Theta(Nlog 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)的二进制翻转。

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

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

当(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=sumlimits_{i=1}^{j-1}dfrac{q_i imes q_j}{(i-j)^2}-sumlimits_{i=j+1}^{n}dfrac{q_i imes q_j}{(i-j)^2},E_i=dfrac{F_i}{q_i})

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

那就化简吧：

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

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

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

[sumlimits_{i=0}^{j}f(i)g(j-i) ]

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

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

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

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

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

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

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

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

两种方法：

令(g'(n-i)=g(i))，原式(=sumlimits_{i=0}^{n-j}f(i+j)g'(n-i)=sumlimits_{i=0}^{n+j}f(i+j)g'(n-i))
令(f'(n-i-j)=f(i+j))，原式(=sumlimits_{i=0}^{n-j}f'(n-i-j)g(i))

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

方法一需要注意的是，为了不影响结果，要保证(sumlimits_{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
", (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
", (q1[i].r - q2[N-i].r)/(double)limit);
	}
	return 0;
}

未完待续...

（五）参考文献

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