FFT快速傅里叶变换详解

介绍

简而言之，这个东西用来做多项式乘法。

比如说，有两个多项式：
$$3x^2+2x+1,2x^2+x+4$$

那么他们乘起来就是
$$6x^4+7x^3+16x^2+9x+4$$

$FFT$ 也就是做了这样一件事。

但是如果暴力取做这个乘法的话，时间复杂度是 $O(n^2)$ 的，然而优秀的 $FFT$ 只需要 $O(nlogn)$。

为了方便，下面设相乘的多项式为 $A,B$，乘起来得到的多项式为 $C$，也就是 $A\times B=C$，以及设 $A$ 是一个 $n$ 次函数，$B$ 是一个 $m$ 次函数，$C$ 是一个 $k$ 次函数。

（有什么不严谨的地方还请见谅qwq）

---

前置芝士

• 三角函数总得要会呀
• 得知道复数是个什么东西
• 平面直角坐标系得会用吧

平面直角坐标系和三角函数就不讲了。要是这个东西都不知道学什么 $FFT$……

复数

前置芝士——虚数： 定义虚数单位 $i=\sqrt{-1}$，虚数就是形如 $a+bi(b\ne 0)$ 的数，其中 $a,b$ 都是实数。

复数： 形如 $a+bi$ 的数，其中 $a,b$ 都是实数。说白了，复数就是实数与虚数的并集。

实部： 也就是 $a$。

虚部： 也就是 $bi$。

复数加法： 实部与实部相加，虚部与虚部相加。 例： $(2+3i)+(-9+4i)=-7+7i$。

复数减法： 与加法类似，实部与实部相减，虚部与虚部相减。

复数乘法：
根据多项式乘法的规则去乘即可。例： $(2-2i)\times (3+1i)=2\times 3+2i-6i-2i^2=6-4i+2=8-4i$。别忘了 $i=\sqrt{-1}$，即 $i^2=-1$。
由此可以得出，复数乘法的规则是：实部与实部相乘 $+$ 虚部与虚部相乘 $=$ 积的实部；实部与虚部相乘 $+$ 虚部与实部相乘 $=$ 积的虚部。

根据这种计算规则，我们可以写出这样的复数结构体：

struct comp{
	double x,y;
	comp(double xx=0,double yy=0):x(xx),y(yy){}
	comp operator +(const comp b){return (comp){x+b.x,y+b.y};}
	comp operator -(const comp b){return (comp){x-b.x,y-b.y};}
	comp operator *(const comp b){return (comp){x*b.x-y*b.y,x*b.y+y*b.x};}
};

---

正题

我们发现，将上面两个多项式乘起来之后，得到了另一个多项式。假如我们将这三个多项式看成三个函数，那么会怎么样呢？

初中时就学过：$n+1$ 个点可以确定一个 $n$ 次函数。那我们只需要找出 $k+1$ 个在函数 $C$ 上面的点，那么就可以求出 $C$ 了。

但是，直接找这些点的过程就是 $O(n^2)$ 的了……

于是，$FFT$ 对这个方面进行了优化！他不找实数点，而是去找复数点。

对于一个复数，我们可以将其表示为 $a+bi$，那么可以在平面直角坐标系中将其表示成一个坐标为 $(a,b)$ 的点。

比如说，对于一个复数 $3+4i$，在平面直角坐标系中就是：

这个点有另外一个表示方法：我们设它与原点之间的线段的长度为 $z$，与 $x$ 轴之间的夹角为 $\theta$，也就是这样子：

因为 $x=\cos (\theta )\times z,y=\sin (\theta )\times z$，所以这个点的坐标又可以表示为：$(\cos (\theta )\times z,\sin (\theta )\times z)$。

单位圆

单位圆： 以平面直角坐标系的原点为圆心做一个半径为 $1$ 的圆，这个圆我们称它为单位圆。

就是这样的一个东西：

单位根

定义： 我们称满足 $x^n=1$ 的复数 $x$ 为 $n$ 次单位根。

当然，对于一个 $n$，是可以有多个单位根的，其中一个就是 $1$。

怎么找呢？

对于一个单位圆，在其内部做一个正 $n$ 边形，要求每个顶点都在圆上，并且有一个顶点在 $(1,0)$ 位置（因为有一个单位根一定是 $1$），那么这个 $n$ 边形的每一个顶点就是一个单位根。

比如说当 $n=6$ 时，图长这个样子：

其中，$w_n$ 表示 $n$ 次单位根，上面标记了 $w_6$ 的点，都是单位根。

我们不妨设 $(0,0)\to (1,0)$ 这条线段逆时针旋转遇到的第一条线段所对应的点为 $w_n$ 的一次方，也就是 $w_n^1$，剩下的点依次为 $w_n^2,w_n^3...$ 也就是：

能这样表示，那么自然也满足幂的运算规则：$x^a\times x^b=x^{a+b}$，在这里也就是 $w_n^a\times w_n^b=w_n^{a+b}$。

证明如下：

我们设原点到 $w_n^1$ 这条线段与 $x$ 轴的夹角为 $\alpha$，发现根据上面的编号规则，其实随着上标的增加，这个夹角也在增加，那么 $w_n^{a+b}$ 的夹角大小就是 $a\times \alpha +b\times \alpha$。

根据上面的坐标的另外一种表示法，可以得到 $w_n^{a+b}$ 的坐标为：$(\cos (a\times \alpha +b\times \alpha )\times 1,\sin (a\times \alpha +b\times \alpha )\times 1)$，也就是 $(\cos (a\times \alpha +b\times \alpha ),\sin (a\times \alpha +b\times \alpha ))$。

考虑两角和公式，有：
$$\begin{gathered} \cos (a\times \alpha +b\times \alpha )=\cos (a\times \alpha )\cos (b\times \alpha )-\sin (a\times \alpha )\sin (b\times \alpha ) \\ \sin (a\times \alpha +b\times \alpha )=\sin (a\times \alpha )\cos (b\times \alpha )+\cos (a\times \alpha )\sin (b\times \alpha ) \end{gathered}$$

这便是 $w_n^{a+b}$ 的坐标。

那么再看 $w_n^a$ 乘 $w_n^b$ 到底乘出了个什么东西。

先写出他们的坐标：
$$\begin{gathered} w_n^a:(\cos (a\times \alpha ),\sin (a\times \alpha )) \\ {w}_n^b:(\cos (b\times \alpha ),\sin (b\times \alpha )) \end{gathered}$$

要记得，这两个坐标相当于他们的实部和虚部。按照复数的运算规则，乘出来的数的坐标跟 $w_n^{a+b}$ 一样，于是得证。

---

由此可以得出，单位根有五大性质：

1. $w_n^0=w_n^n=1$
2. $w_n^0,w_n^1,...,w_n^{n-1}$ 互不相同
3. $w_n^k=w_{2n}^{2k}$
4. $w_n^k=-w_n^{k+\frac{n}{2}}(k<\frac{n}{2})$
5. ${\sum }_{i=0}^{n-1}(w_n^{j-k}{)}^i=\{\begin{array}{l}0(j\ne k)\\ n(j=k)\end{array}$

性质1证明： 过于明显就不证了。

性质2证明： 过于明显依然不证。

性质3证明：

感性理解：下标从 $n$ 变成了 $2n$，意味着单位圆被分成了原来的两倍份，这样的话每一份的大小只有原来的 $\frac{1}{2}$，那么因为我取的份数变成了原来的两倍，两倍乘 $\frac{1}{2}$ 等于 $1$，也就是和原来一样。

理性证明：带入定义式即可。
$$\begin{gathered} w_{2n}^{2k}=\cos (\frac{360°}{2n}\times 2k)+\sin (\frac{360°}{2n}\times 2k)i \\ =\cos (\frac{360°}{n}\times k)+\sin (\frac{360°}{n}\times k)i \\ =w_n^k \end{gathered}$$

性质4证明：

感性理解：转半圈就变成原来的相反数了。

理性证明：带入定义式即可。
$$\begin{gathered} -w_n^{k+\frac{n}{2}}=-w_n^k\times w_n^{\frac{n}{2}} \\ =-w_n^k\times (\cos (\frac{360°}{n}\times \frac{n}{2})+\sin (\frac{360°}{n}\times \frac{n}{2})i) \\ =-w_n^k\times (\cos (180°)+\sin (180°)i) \\ =-w_n^k\times (-1+0) \\ =w_n^k \end{gathered}$$

性质5证明： 难以感性理解，直接上证明。

当 $j\ne k$ 时，考虑等比数列求和。
$$\begin{gathered} \sum _{i=0}^{n-1}(w_n^{j-k}{)}^i=\frac{1-(w_n^{j-k}{)}^n}{1-w_n^{j-k}} \\ =\frac{1-(w_n^n{)}^{j-k}}{1-w_n^{j-k}} \\ =\frac{1-1}{1-w_n^{j-k}} \\ =0 \end{gathered}$$

当 $j=k$ 时，有：
$$\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i=\sum _{i=0}^{n-1}(w_n^0{)}^i=\sum _{i=0}^{n-1}1=n$$

于是得证。

务必记得这些性质，当你看不懂下面的柿子时，就回来看看吧。

---

弧度制

在 $C++$ 里面，你是无法输入一个度数的，于是你需要弧度制。

具体来说，弧度制就是用一个弧来表示一个度数。

对于每一个度数 $x°$，他有一个对应的弧长，这段弧长的长度就是单位圆中做一个度数为 $x°$ 的圆心角，这个圆心角所对弧的长度。

打个比方，比如说 $60°$ 对应的弧长。

单位圆的半径是 $1$，那么周长就是 $2\pi$，那么这段弧的长度就是 $\frac{2\pi }{6}=\frac{1}{3}\pi$。

在 $C++$ 中，$\frac{1}{3}\pi$ 就相当于 $60°$。

那么怎么在 $C++$ 中求 $\pi$ 呢？

我们注意到，单位圆周长的一半就是 $\pi$，也就是说，$\pi$ 对应 $180°$，因为 $\cos (180°)=-1$，所以，我们可以用 $acos(-1)$ 来算出 $\pi$。（注：$acos$ 也就是 $\cos$ 的逆运算，$C++$ 中是有这个函数的）

---

快速傅里叶变换

回到上面的话题，$FFT$ 去找复数点，并不是指找复数平面直角坐标系上的点，而是去找以复数作为横坐标的点，因为复数平面直角坐标系上的点表示的终究只是一个数，而我们的最终目的是要找若干个真正的点。

要知道，复数既然是一个数，那么将这个数带入到函数式里面，也肯定是能得到一个值的，现在我们要做的，就是求出将 $w_n^0,w_n^1,...,w_n^{n-1}$ 带入函数 $C$ 中的值。

但是我们并不知道函数 $C$ 的解析式呀，怎么求呢？

其实我们只需要求出将 $x$ 带入函数 $A$ 中的值 $y_1$ 以及带入函数 $B$ 中的值 $y_2$，然后乘起来就能得到将 $x$ 带入函数 $C$ 中的值，也就是 $y_1\times y_2$。

现在问题变成了如何快速求出函数 $A$ 和函数 $B$ 在 $w_n^0,w_n^1,...w_n^{n-1}$ 处的取值。

因为这两个问题是一样的，我们下面只考虑如何求出函数 $A$ 在 $w_n^0,w_n^1,...w_n^{n-1}$ 处的取值。

我们称，函数在 $w_n^0,w_n^1,...w_n^{n-1}$ 处的取值为这个函数的傅里叶变换。

我们设函数 $A$ 的解析式为 ${\sum }_{i=0}^{n-1}{a}_i{x}^i$，将这 $n$ 项根据 $x$ 的上标的奇偶性分类，变成：
$$A(x)=(a_0+a_2{x}^2+a_4{x}^4+...)+(a_{1}x+a_3{x}^3+a_5{x}^5+...)$$

设 $A_1(x)=a_0+a_{2}x+a_4{x}^2+...,A_2(x)=a_1+a_{3}x+a_5{x}^2+...$，那么有：
$$A(x)=A_1(x^2)+x{A}_2(x^2)$$

我们考虑将 $n$ 个单位根分成两部分来做，一部分是 $w_n^k(k<\frac{n}{2})$，另一部分是 $w_n^{k+\frac{n}{2}}(k<\frac{n}{2})$。

先将 $w_n^k$ 带入：
$$\begin{gathered} A(w_n^k)=A_1((w_n^k{)}^2)+w_n^k{A}_2((w_n^k{)}^2) \\ =A_1(w_{\frac{n}{2}}^k)+w_n^k{A}_x(w_{\frac{n}{2}}^k) \end{gathered}$$

再将 $w_n^{k+\frac{n}{2}}$ 带入：
$$\begin{gathered} A(w_n^{k+\frac{n}{2}})=A_1((w_n^{k+\frac{n}{2}}{)}^2)+w_n^{k+\frac{n}{2}}{A}_x((w_n^{k+\frac{n}{2}}{)}^2) \\ =A_1((-w_n^k{)}^2)-w_n^k{A}_x((-w_n^k{)}^2) \\ =A_1(w_n^{2k})-w_n^k{A}_x(w_n^{2k}) \\ =A_1(w_{\frac{n}{2}}^k)-w_n^k{A}_x(w_{\frac{n}{2}}^k) \end{gathered}$$

诶？上下两个柿子只差一个符号！

这说明，当我们求出 $A_1$ 和 $A_2$ 时，就可以求出 $A$ 了。

于是问题从 求 $A$ 在 $w_n^0,w_n^1,...,w_n^{n-1}$ 处的取值 变成了求 $A_1$ 和 $A_2$ 在 $w_{\frac{n}{2}}^0,w_{\frac{n}{2}}^1,...,w_{\frac{n}{2}}^{\frac{n}{2}-1}$ 处的取值。问题的规模缩小了一半，这样就可以愉快地分治了！

于是代码如下：

void fft(comp *f,int len)
{
	if(len==1)return;//假如只有常数项，直接返回
	comp A1[len/2+1],A2[len/2+1];
	
	for(int i=0;i<(len>>1);i++)//按照奇偶性分类
	A1[i]=f[i*2],A2[i]=f[i*2+1];
	
	fft(A1,len>>1);fft(A2,len>>1);//递归下去
	
	comp wn(cos(2*pi/len),sin(2*pi/len)),w(1,0);
	//求出wn1(即wn)，w用来做累乘，求出wnk
	for(int i=0;i<(len>>1);i++,w=w*wn)
	f[i]=A1[i]+A2[i]*w,f[i+(len>>1)]=A1[i]-A2[i]*w;//求解
}

做完了？不不不，滚动条告诉你并没有这么简单。

我们现在只是求出了这些点，我们还没得出解析式呢。

不对，点都还没求出来，还差这么一步：

fft(a,up);fft(b,up);
for(int i=0;i<up;i++)
c[i]=a[i]*b[i];

我们现在求出的 $c$ 数组，相当于函数 $C$ 的傅里叶变换，怎么通过这个傅里叶变换得到函数 $C$呢？

快速傅里叶逆变换

其实和快速傅里叶变换差不多，快速傅里叶变换是函数在 $w_n^0,w_n^1,...w_n^{n-1}$ 处的取值，而逆变换就是在 $w_n^0,w_n^{-1},...w_n^{-n+1}$ 处的取值。

显然 $w_n^{-1}=w_n^{n-1}$，因为 $w_n^{-1}=w_n^{-1}\times 1=w_n^{-1}\times w_n^n=w_n^{n-1}$。

这个很好实现，只需要将上面代码里的 $w_n^1$ 换成 $w_n^{-1}$ 即可，相当于顺时针旋转。发现 $w_n^1$ 和 $w_n^{-1}$ 的横坐标相同，纵坐标互为相反数，所以只需要将上面的

wn(cos(2*pi/len),sin(2*pi/len))

改成

wn(cos(2*pi/len),-sin(2*pi/len))

回到正题

我们设 $d$ 为 $c$（这里的 $c$ 是上面求出的 $c$ 数组，也就是函数 $C$ 的傅里叶变换）的傅里叶逆变换，那么有：
$$\begin{gathered} d_k=\sum _{i=0}^{n-1}{c}_i(w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{C}_j(w_n^i{)}^j)(w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{C}_j(w_n^j{)}^i(w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{C}_j(w_n^j\times w_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{C}_j(w_n^{j-k}{)}^i \\ =\sum _{j=0}^{n-1}{C}_j\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i \end{gathered}$$

这时候我们需要掏出我们封尘已久的性质5！

性质5：${\sum }_{i=0}^{n-1}(w_n^{j-k}{)}^i=\{\begin{array}{l}0(j\ne k)\\ n(j=k)\end{array}$

带入得：
$$\begin{gathered} =\sum _{j=0}^{n-1}{C}_j\times \{\begin{array}{l}0(j\ne k)\\ n(j=k)\end{array} \\ =C_k\times n \end{gathered}$$

诶嘿？

也就是说，我们求出 $c$ 数组的快速傅里叶逆变换后，将每一位除以 $n$，就可以得到 $C$ 函数了！

你难道又以为完了？不不不，还有个细节问题。

发现我们每次需要将当前的序列进行奇偶性分类，分成两份，然后再合并得到当前序列的答案，那么这就要求，当前序列的长度是偶数。进而要求，原序列的长度需要是 $2$ 的某个幂。

所以我们需要找到第一个大于序列长度的 $2$ 的幂，让它来做序列长度，比原来多的位置不用管，放 $0$ 即可。

代码是这样的：

while(up<=n+m)up<<=1;//n是函数A的长度，m是函数B的长度
//因为我们需要让up大于函数C的长度，所以需要让up>n+m

然后先给出逆变换的代码：

void fft(comp *f,int len,int type)//type表示做的是变换还是逆变换
{
	if(len==1)return;
	comp A1[len/2+1],A2[len/2+1];
	for(int i=0;i<(len>>1);i++)
	A1[i]=f[i*2],A2[i]=f[i*2+1];
	fft(A1,len>>1,type);fft(A2,len>>1,type);
	comp wn(cos(2*pi/len),type*sin(2*pi/len)),w(1,0);//注意这里type要乘上纵坐标
	for(int i=0;i<(len>>1);i++,w=w*wn)
	f[i]=A1[i]+A2[i]*w,f[i+(len>>1)]=A1[i]-A2[i]*w;
}

模板题 完整代码：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 3000010

int n,m;
struct comp{
	double x,y;
	comp(double xx=0,double yy=0):x(xx),y(yy){}
	comp operator +(const comp b){return (comp){x+b.x,y+b.y};}
	comp operator -(const comp b){return (comp){x-b.x,y-b.y};}
	comp operator *(const comp b){return (comp){x*b.x-y*b.y,x*b.y+y*b.x};}
};
comp a[maxn],b[maxn],c[maxn];
int up=1;
const double pi=acos(-1);
void fft(comp *f,int len,int type)
{
	if(len==1)return;
	comp A1[len/2+1],A2[len/2+1];
	for(int i=0;i<(len>>1);i++)
	A1[i]=f[i*2],A2[i]=f[i*2+1];
	fft(A1,len>>1,type);fft(A2,len>>1,type);
	comp wn(cos(2*pi/len),type*sin(2*pi/len)),w(1,0);
	for(int i=0;i<(len>>1);i++,w=w*wn)
	f[i]=A1[i]+A2[i]*w,f[i+(len>>1)]=A1[i]-A2[i]*w;
}

int main()
{
	scanf("%d %d",&n,&m);
	for(int i=0;i<=n;i++)
	scanf("%lf",&a[i].x);
	for(int i=0;i<=m;i++)
	scanf("%lf",&b[i].x);
	while(up<=n+m)up<<=1;
	fft(a,up,1);fft(b,up,1);
	for(int i=0;i<up;i++)
	c[i]=a[i]*b[i];
	fft(c,up,-1);
	for(int i=0;i<=n+m;i++)
	printf("%d ",(int)(c[i].x/(double)up+0.5));
}

很抱歉的告诉你，$FFT$ 还没完。

这份优秀的代码提交之后，你可以获得 $80$ 分的好成绩。

为啥？

发现在分治的每一层我们都开了两个新的 $A_1$ 和 $A_2$ 数组，也就是说，我们有 $nlogn$ 的空间是开在栈里面的，这就直接导致了栈空间爆炸。

于是，当数据大起来时，除了 $AC$ 之外，你还可以获得 $TLE$ 或 $RE$。

于是，有了下一节。

迭代FFT

我们不能每次都把它进行奇偶性分类，那么不如一次分到底，然后从下往上合并？

假设我们已经分完了，那么代码就是这样的了：

void fft(comp *f,int len,int type)
{
	for(int mid=1;mid<len;mid<<=1)//枚举一块的大小
	{
		comp wn(cos(pi/mid),type*sin(pi/mid));//求出单位根
		for(int block=mid<<1,j=0;j<len;j+=block)
		//block是两块的大小，因为我们需要将相邻两块合成一整块，j枚举所有整块
		{
			comp w(1,0);
			for(int i=j;i<j+mid;i++,w=w*wn)
			//i枚举这一整块中的左边那一块，i+mid也就是在枚举右边那一块
			{
				comp x=f[i],y=f[i+mid]*w;//合并两块
				f[i]=x+y;f[i+mid]=x-y;
			}
		}
	}
}

是不是更简单啦~

剩下只有一个问题，如何求出每个数最后会被分到的位置？

找规律大法好！

观察一下每个数和他最后的位置的二进制数的关系。（第一行是这个数的二进制数，第二行是这个位置的二进制数）

相信大家肯定动动脚指头就能发现，每个数的二进制数的翻转就是他最后的位置。

那么暴力求即可。

---

我猜你又以为完了。

事实上并不是。

这个东西还有 $O(n)$ 的求法~

对于一个数 $i$，我们将他的最后一位拿走，然后进行翻转，再把拿走的那一位放在第一位，就可以得到 $i$ 的翻转。这就是 $O(n)$ 的求法，因为把他的最后一位拿走之后的翻转我们已经知道了，所以可以直接得到答案。

打个比方：对于一个数 $181$，他的二进制数是 $10110101$。

那么先拿走最后一位，变成：
$$1011010$$

翻转，变成
$$0101101$$

将一开始拿走的 $1$ 放到第一位变成：
$$10101101$$

于是我们得到了 $i$ 的翻转。

怎么知道 $i$ 被拿走最后一位的那个数的翻转呢？

$i$ 被拿走最后一位，相当于 $\frac{i}{2}$，那么这个翻转就是 $\frac{i}{2}$ 的翻转除以 $2$。

为什么要除以 $2$？

因为 $\frac{i}{2}$ 在翻转的时候，最高位被补了 $0$，翻转之后这个 $0$ 就跑到了最低位，所以我们要除以 $2$ 把它去掉。

依然用上面的那个栗子：

拿走最后一位之后，$i$ 变成
$$1011010$$

最高位补 $0$，变成：
$$01011010$$

翻转，变成：
$$01011010$$

除以 $2$，变成：
$$0101101$$

这才是我们想要的。

那么问题来了，为什么 $\frac{i}{2}$ 的最高位会被补 $0$ ？

我们发现，在那个图片里，$001$ 翻转变成 $100$，而不是 $1$ 翻转变成 $1$，这说明，所有数被翻转的位数是固定的，不是根据这个数本身的位数而定的。而这个固定的位数，就是最大的那个数的位数。

于是我们可以得到递推式：

for(int i=1;i<up;i++)//第0位做了跟不做没区别
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
//r[i>>1]>>1 也就是去掉最低位后的翻转
//(i&1) 得到最低位，左移(l-1)位之后就是最高位了，然后与(r[i>>1]>>1)合并

模板题 $AC$ 代码：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 4000010

const double pi=acos(-1);
int n,m;
struct comp{
	double x,y;
	comp(double xx=0,double yy=0):x(xx),y(yy){}
	comp operator +(const comp b){return (comp){x+b.x,y+b.y};}
	comp operator -(const comp b){return (comp){x-b.x,y-b.y};}
	comp operator *(const comp b){return (comp){x*b.x-y*b.y,x*b.y+y*b.x};}
};
comp a[maxn],b[maxn];
int r[maxn];
void fft(comp *f,int len,int type)
{
	for(int i=1;i<len;i++)//这一步下面再解释
	if(i<r[i])swap(f[i],f[r[i]]);
	for(int mid=1;mid<len;mid<<=1)
	{
		comp w(cos(pi/mid),type*sin(pi/mid));
		for(int block=mid<<1,j=0;j<len;j+=block)
		{
			comp p(1,0);
			for(int i=j;i<j+mid;i++,p=p*w)
			{
				comp x=f[i],y=f[i+mid]*p;
				f[i]=x+y;f[i+mid]=x-y;
			}
		}
	}
}
void FFT(comp *p1,comp *p2,int len)
{
	fft(p1,len,1);fft(p2,len,1);
	for(int i=0;i<len;i++)
	p1[i]=p1[i]*p2[i];
	fft(p1,len,-1);
}

int main()
{
	scanf("%d %d",&n,&m);
	for(int i=0;i<=n;i++)
	scanf("%lf",&a[i].x);
	for(int i=0;i<=m;i++)
	scanf("%lf",&b[i].x);
	int up=1,l=0;
	while(up<=n+m)up<<=1,l++;
	for(int i=1;i<up;i++)
	r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
	FFT(a,b,up);
	for(int i=0;i<=n+m;i++)
	printf("%d ",(int)(a[i].x/up+0.5));
}

发现代码里面有一个好像很奇怪的地方：

if(i<r[i])swap(f[i],f[r[i]]);

为什么需要判断 $i<r[i]$ 呢？

因为 $r[i]$ 是 $i$ 的翻转，假如有 $r[a]=b$，那么肯定有 $r[b]=a$。

假如我不加那个判断，就变成了：

swap(f[i],f[r[i]]);

那么当你枚举到 $a$ 时，你会将 $a,b$ 交换，枚举到 $b$ 时，你又会将 $b,a$ 交换，那么不就相当于没有变化？

所以，不妨设 $a<b$，那么加上那个判断之后，只有在枚举到 $a$ 的时候会进行交换，而枚举到 $b$ 时就不会交换了。

恭喜你，学完 $FFT$ 了！

一点题表

ZJOI 2014 力   题解
AH2017/HNOI2017 礼物   题解