FFT-绝对能看懂的-傅里叶变换-掰开揉碎解释-三个c++代码

tip:本文由typora复制过来并修改了一些显示错误，如果还有显示不正确，或者文章内容有误的地方，欢迎在评论区指出

快速傅里叶变换(FFT)详解

前言

DFT：离散傅里叶变换—>O(n^2)计算多项式乘法

FFT：快速傅里叶变换—>O(n∗log⁡(n))计算多项式乘法

FNTT/NTT：快速傅里叶变换的优化版—>优化常数及误差

FWT：快速沃尔什变换—>利用类似FFT的东西解决一类卷积问题

MTT：毛爷爷的FFT—>非常nb/任意模数

FMT 快速莫比乌斯变化

多项式

系数表示法

设A(x)表示一个n−1次多项式

则
$$A(x)=\sum _{i=1}^n{a}_i\ast x^i$$
利用这种方法计算多项式乘法复杂度为O(n)（第一个多项式中每个系数都需要与第二个多项式的每个系数相乘）
$$例如：A(3)=2+3\ast x+x^2$$

点值表示法

将n互不相同的x带入多项式，会得到n个不同的取值y

则该多项式被这n个点$(x_1,y_1),(x_2,y_2),\dots ,(x_n,y_n)$唯一确定

其中
$$y_i=\sum _{j=0}^{n-1}{a}_j\ast x_i^j$$
例如：上面的例子用点值表示法可以为(0,2),(1,5),(2,12)

利用这种方法计算多项式乘法的时间复杂度仍然为
$$O(n^2)（选点O(n)，每次计算O(n)）$$

---

我们可以看到，两种方法的时间复杂度都为O*(*n^2)，我们考虑对其进行优化

对于第一种方法，由于每个点的系数都是固定的，想要优化比较困难

对于第二种方法，貌似也没有什么好的优化方法，不过当你看完下面的知识，或许就不这么想了

---

复数

在介绍复数之前，首先介绍一些可能会用到的东西

向量

同时具有大小和方向的量

在几何中通常用带有箭头的线段表示

圆的弧度制

等于半径长的圆弧所对的圆心角叫做1弧度的角，用符号rad表示，读作弧度。用弧度作单位来度量角的制度叫做弧度制

公式:
$$\begin{gathered} 1^{\circ }=\frac{\pi }{180} \\ 180^{\circ }=\pi \ast rad \end{gathered}$$

平行四边形定则

平行四边形定则：AB+AD=AC

复数

定义

设a,b为实数，i²=-1，形如a+bi的数叫负数，其中ii被称为虚数单位，复数域是目前已知最大的域

在复平面中，x代表实数，y轴（除原点外的点）代表虚数，从原点(0,0)到(a,b)的向量表示复数a+bi

模长：从原点(0,0)到点(a,b)的距离，即$\sqrt{a^2+b^2}$

幅角：假设以逆时针为正方向，从xx轴正半轴到已知向量的转角的有向角叫做幅角

运算法则

加法：

因为在复平面中，复数可以被表示为向量，因此复数的加法与向量的加法相同，都满足平行四边形定则（就是上面那个）

乘法：

几何定义：复数相乘，模长相乘，幅角相加

代数定义：
$$\begin{gathered} (a+bi)\ast (c+di) \\ =ac+adi+bci+bd{i}^2 \\ =ac+adi+bci-bd \\ =(ac-bd)+(bc+ad)i \end{gathered}$$

单位根

下文中，默认$n$为$2$的正整数次幂
在复平面上，以原点为圆心，1为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的$n$等分点为终点，
做n个向量，设幅角为正且最小的向量对应的复数为${\omega }_n$，称为$n$次单位根。
根据复数乘法的运算法则，其余n-1个复数为${\omega }_n^2,{\omega }_n^3,\dots ,{\omega }_n^n$
注意${\omega }_n^0={\omega }_n^n=1$（对应复平面上以x轴为正方向的向量）

那么如何计算它们的值呢？这个问题可以由欧拉公式解决
$${\omega }_n^k=\cos k\ast \frac{2\pi }{n}+i\sin k\ast \frac{2\pi }{n}$$
单位根的幅角为周角的$\frac{1}{n}$

例如

图中向量AB表示的复数为8次单位根

• 在代数中，若$z^n=1$，我们把z称为n次单位根

单位根的性质

$1.{\omega }_n^k=cosk\frac{2\pi }{n}（即上面的公式）$
证明：这个第一步到第二步由定义得出，第二步到第三步由欧拉公式得出

$2.{\omega }_{2n}^{2k}={\omega }_n^k$
证明：${\omega }_{2k}^{2n}=e^{\frac{2\pi 2ki}{2n}}=e^{\frac{2\pi ki}{n}}={\omega }_n^k$

$3.{\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
证明：$$\begin{gathered} {\omega }_n^{k+\frac{n}{2}}={\omega }_{2n}^{2k}\ast {\omega }_{2n}^n \\ ={\omega }_n^k\ast (cos\pi +i\ast sin\pi )=-{\omega }_n^k \end{gathered}$$

$4.\omega 0n=\omega nn=1$
证明：不用证了吧……

---

正文之前

这段话有可能有助于您理解本算法：

傅立叶这个大神仙根本就没见过计算机长什么样，所以他提出的傅立叶变换和逆变换只是一种将系数转点值和将点值转系数的方法，没有任何降低复杂度的功效，至于快速傅立叶变换是后人再研究傅立叶变换发现的一种加速方法，是对$DFT$和$IDFT$的优化

离散傅立叶变换（DFT）

假设$f(x)={\sum }_{i=0}^{n-1}{a}_i\ast x^i$

$DFT(a)=(f(1),f({\omega }_n),f({\omega }_n^2),\dots \dots ,f({\omega }_n^{n-1}))$

通俗点说，就是对于一个系数表示法的多项式$f(x)$，将$(1,{\omega }_n,{\omega }_n^2,\dots \dots ,{\omega }_n^{n-1})$带入求出该多项式的点值表示法

离散傅立叶逆变换（IDFT）

将$f(x)$在$n$个$n$次单位根处的点值表示转化为系数表示

这里就可以回答，为什么我们要让$n$次单位根作为$x$代入多项式

假设$(y_0,y_1,y_2,\dots \dots ,y_{n-1})$是多项式$A(x)={\sum }_{i=1}^n{a}_i\ast x^i$的离散傅立叶变换。

我们另有一个多项式$B(x)={\sum }_{i=1}^n{y}_i\ast x^i$

将上述$n$次单位根的倒数$(1,{\omega }^{-}{1}_n,{\omega }_n^{-2},\dots \dots ,{\omega }_n^{-(n-1)})$带入$B(x)$

得到新的离散傅立叶变换$(z_0,z_1,z_2,\dots \dots ,z_{n-1})$

则我们发现
$$\begin{gathered} z_k=\sum _{i=0}^{n-1}{y}_i\ast ({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j\ast ({\omega }_n^i{)}^j)\ast ({\omega }_n^{-k}{)}^i \\ =\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i) \end{gathered}$$
我们单独考虑：

当$j-k=0$时 , 答案为$n$

当$j\ne k$时 , 等比数列求和得到

$\frac{({\omega }_n^{j-k}{)}^n-1}{{\omega }_n^{j-k}-1}=\frac{({\omega }_n^n{)}^{j-k}-1}{{\omega }_n^{j-k}-1}=\frac{1^{j-k}-1}{{\omega }_n^{j-k}-1}=0$

所以

$\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)=n\ast a_k$

即

$a_k=\frac{z_k}{n}$

得出结论：对于以$A(x)$的离散傅立叶变换作为系数的多项式$B(x)$，取单位根的倒数$(1,{\omega }_n^{-1},{\omega }_n^{-2},\dots \dots ,{\omega }_n^{-(n-1)})$作为$x$代入，再将结果除以$n$即为$A(x)$的系数

这个结论实现了将多项式点值转化为系数

---

讲了这么多，貌似跟我们的正题没啥关系啊。。

OK！各位坐稳了，前方高能！

---

快速傅里叶变换FFT

我们前面提到过，一个n次多项式可以被n个点唯一确定。

那么我们可以把单位根的0到n−1次幂带入，这样也可以把这个多项式确定出来。但是这样仍然是O(n^2)的呀！

我们设多项式A(x)的系数为(a0,a1,a2,…,an−1),即
$$A(x)=a_0+a_{1}x+a_2{x}^2+a_3\ast x^3+a_4\ast x^4+a_5\ast x^5+\cdots +a_{n-2}\ast x^{n-2}+a_{n-1}\ast x^{n-1}$$
将其下标按照奇偶性分类
$$A(x)=(a_0+a_2\ast x^2+a_4\ast x^4+\cdots +a_{n-2}\ast x^{n-2})+(a_1\ast x+a_3\ast x^3+a_5\ast x^5+\cdots +a_{n-1}\ast x^{n-1})$$
我们设
$$\begin{gathered} A_1(x)=a_0+a_2\ast x+a_4\ast x^2+\cdots +a_{n-2}\ast x^{\frac{n}{2}-1} \\ {A}_2(x)=a_1\ast x+a_3\ast x+a_5\ast x^2+\cdots +a_{n-1}\ast x^{\frac{n}{2}-1} \end{gathered}$$
那么不难得到
$$A(x)=A_1(x^2)+x\ast A_2(x^2)$$

$$我们将x={\omega }_n^k(k<\frac{n}{2})代入得$$

$$\begin{gathered} A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})式子① \\ =A_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$

$$同理，将x={\omega }_n^{k+\frac{n}{2}}代入得$$

$$\begin{array}{rl}& A({\omega }_n^{k+\frac{n}{2}})=A_1({\omega }_n^{2k}\ast {\omega }_n^n)-{\omega }_n^k{A}_2({\omega }_n^{2k}\ast {\omega }_n^n)\\ & =A_1({\omega }_n^{2k})-{\omega }_n^k{A}_2({\omega }_n^{2k})式子②\\ & =A_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k)\end{array}$$

大家有没有发现什么规律？

没错！这两个式子只有一个常数项不同！

那么当我们在枚举第一个式子的时候，我们可以O(1)的得到第二个式子的值

又因为第一个式子的k在取遍[0,n/2−1]时，k+n/2取遍了[n/2,n−1]

所以我们将原来的问题缩小了一半！

$$\begin{gathered} 只要求出A_1(x)和A_2(x)在({\omega }_{\frac{n}{2}}^0,{\omega }_{\frac{n}{2}}^1,\dots \dots ,{\omega }_{\frac{n}{2}}^{\frac{n}{2}-1})的点值表示， \\ 就可以O(n)地求出A(x)在(1,{\omega }_n,{\omega }_n^2,\dots \dots ,{\omega }_n^{n-1}) \end{gathered}$$

而缩小后的问题仍然满足原问题的性质，所以我们可以递归的去搞这件事情！

直到多项式仅剩一个常数项，这时候我们直接返回就好啦

时间复杂度：

不难看出FFT是类似于线段树一样的分治算法。

因此它的时间复杂度为$O(nlogn)$

快速傅里叶逆变换

不要以为FFT到这里就结束了。

我们上面的讨论是基于点值表示法的。

但是在平常的学习和研究中很少用点值表示法来表示一个多项式。

所以我们要考虑如何把点值表示法转换为系数表示法，这个过程叫做傅里叶逆变换

$我们假设(y_0,y_1,y_2,\dots ,y_{n-1})为(a_0,a_1,a_2,\dots ,a_{n-1})的傅里叶变换（即点值表示）$

那么变换向量$(c_0,c_1,c_2,\dots ,c_{n-1})$满足
$$c_k=\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i$$

即$多项式B(x)=y_0,y_{1}x,y_2{x}^2,\dots ,y_{n-1}{x}^{n-1}$在${\omega }_n^0,{\omega }_n^{-1},{\omega }_n^{-2},\dots ,{\omega }_{n-1}^{-(n-1)}$处的点值表示

$(c_0,c_1,c_2,\dots ,c_{n-1})满足$
$$\begin{gathered} c_k=\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j)({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i)({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i) \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i \\ =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i \\ =\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i) \end{gathered}$$
设$S(x)={\sum }_{i=0}^{n-1}{x}^i$
将${\omega }_n^k$代入得
$S({\omega }_n^k)=1+({\omega }_n^k)+({\omega }_n^k{)}^2+\dots ({\omega }_n^k{)}^{n-1}$
当k!=0时

等式两边同乘${\omega }_n^k$得
$${\omega }_n^{k}S({\omega }_n^k)={\omega }_n^k+({\omega }_n^k{)}^2+({\omega }_n^k{)}^3+\cdots +({\omega }_n^k{)}^n$$
两式相减得
$$\begin{gathered} {\omega }_n^{k}S({\omega }_n^k)-S({\omega }_n^k)=({\omega }_n^k{)}^n-1 \\ S({\omega }_n^k)=\frac{({\omega }_n^k{)}^n-1}{{\omega }_n^k-1} \\ S({\omega }_n^k)=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1} \\ S({\omega }_n^k)=\frac{1-1}{{\omega }_n^k-1} \end{gathered}$$

观察这个式子，不难看出它分母不为0，但是分子为0

因此，当K!=0时，$S({\omega }_n^k)=0$

那当k=0时呢？

很显然，
$$S({\omega }_n^0)=n$$
继续考虑刚刚的式子
$$c_k=\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)$$
当$j\ne k时，值为0$
当$j=k时，值为n$
因此，
$$
c_k=na_k\

a_k=\frac{c_k}{n}
$$
这样我们就得到点值与系数之间的表示啦

理论总结

至此，FFT的基础理论部分就结束了。

我们来小结一下FFT是怎么成功实现的

首先，人们在用系数表示法研究多项式的时候遇阻

于是开始考虑能否用点值表示法优化这个东西。

然后根据复数的两条性质（这个思维跨度比较大）得到了一种分治算法。

最后又推了一波公式，找到了点值表示法与系数表示法之间转换关系。

emmmm

其实FFT的实现思路大概就是

系数表示法—>点值表示法—>系数表示法

引用一下远航之曲大佬的图

递归版 FFT

递归实现的方法比较简单。

就是按找我们上面说的过程，不断把要求的序列分成两部分，再进行合并

在c++的STL中提供了现成的complex类，但是我不建议大家用，毕竟手写也就那么几行，而且万一某个毒瘤卡STL那岂不是很ＧＧ？

#include<bits/stdc++.h>
using namespace std;
#define mid ((l+r)>>1)
#define eps (1e-8)
inline int read()
{
	int x = 0;
	char ch, f = 1;
	for (ch = getchar(); (ch < '0' || ch > '9') && ch != '-'; ch = getchar());
	if (ch == '-') f = 0, ch = getchar();
	while (ch >= '0' && ch <= '9')
	{
		x = (x << 1) + (x << 3) + ch - '0';
		ch = getchar();
	}
	return f ? x : -x;
}
const int N = 5e6 + 10;
const double pi = acos(-1.0);
int n, m;
struct Complex
{
	double x, y;
	Complex(double tx = 0, double ty = 0)
	{
		x = tx, y = ty;
	}
	inline Complex operator + (const Complex t) const
	{
		return Complex(x + t.x, y + t.y);
	}
	inline Complex operator - (const Complex t) const
	{
		return Complex(x - t.x, y - t.y);
	}
	inline Complex operator * (const Complex t) const
	{
		return Complex(x * t.x - y * t.y, x * t.y + y * t.x);//不懂的看复数的运算那部分
	}
} a[N], b[N];
inline void fft(int limit, Complex *a, int inv)
{
	if (limit == 1) return;//只有一个常数项
	Complex a1[limit >> 1], a2[limit >> 1];
	for (int i = 0; i < limit; i += 2)//根据下标的奇偶性分类
	{
		a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
	}
	fft(limit >> 1, a1, inv);
	fft(limit >> 1, a2, inv);
	Complex Wn = Complex(cos(2.0 * pi / limit), inv * sin(2.0 * pi / limit)), w = Complex(1, 0);//Wn为单位根，w表示幂
	for (int i = 0; i < (limit >> 1); ++i, w = w * Wn) //这里的w相当于公式中的k
	{
		a[i] = a1[i] + w * a2[i];
		a[i + (limit >> 1)] = a1[i] - w * a2[i];//利用单位根的性质，O(1)得到另一部分
	}
}
signed main()
{
	n = read(), m = read();
	for (int i = 0; i <= n; ++i) a[i].x = read();
	for (int i = 0; i <= m; ++i) b[i].x = read();
	int limit = 1;
	while (limit <= n + m) limit <<= 1;
	fft(limit, a, 1);
	fft(limit, b, 1);
    //后面的1表示要进行的变换是什么类型
    //1表示从系数变为点值
    //-1表示从点值变为系数
    //至于为什么这样是对的，可以参考一下c向量的推导过程，
	for (int i = 0; i <= limit; ++i)
	{
		a[i] = a[i] * b[i];
	}
	fft(limit, a, -1);
	for (int i = 0; i <= n + m; ++i) printf("%d ", (int)(a[i].x / limit + 0.5)); //按照我们推倒的公式，这里还要除以n
	return 0;
}

然而我们发现好像有点慢

迭代优化

众所周知递归比较慢，我们有没有什么方法可以用迭代代替递归呢？

通过一顿找规律，我们根本不能发现，每个数字在分治后的位置就是它所在位置的二进制翻转

这个规律也有一个好听的名字，叫蝴蝶定理

那么我们只要预处理出每个数字在最后一次递归中的位置，然后自底向上合并，岂不是可以摆脱递归

我们需要求的序列实际是原序列下标的二进制反转！

因此我们对序列按照下标的奇偶性分类的过程其实是没有必要的

这样我们可以O(n)的利用某种操作得到我们要求的序列，然后不断向上合并就好了

#include<bits/stdc++.h>
using namespace std;
#define eps (1e-8)
inline int read()
{
	int x = 0;
	char ch, f = 1;
	for (ch = getchar(); (ch < '0' || ch > '9') && ch != '-'; ch = getchar());
	if (ch == '-') f = 0, ch = getchar();
	while (ch >= '0' && ch <= '9')
	{
		x = (x << 1) + (x << 3) + ch - '0';
		ch = getchar();
	}
	return f ? x : -x;
}
const int N = 5e6 + 10;
const double pi = acos(-1.0);
int n, m;
int limit = 1, len;
int pos[N];
struct Complex
{
	double x, y;
	Complex(double tx = 0, double ty = 0)
	{
		x = tx, y = ty;
	}
	inline Complex operator + (const Complex t) const
	{
		return Complex(x + t.x, y + t.y);
	}
	inline Complex operator - (const Complex t) const
	{
		return Complex(x - t.x, y - t.y);
	}
	inline Complex operator * (const Complex t) const
	{
		return Complex(x * t.x - y * t.y, x * t.y + y * t.x);
	}
} a[N], b[N], buf[N];
inline void fft(Complex *a, int inv)
{
	for (int i = 0; i < limit; ++i)
		if (i < pos[i]) swap(a[i], a[pos[i]]); //求出要迭代的序列
	for (int mid = 1; mid < limit; mid <<= 1) //待合并区间的长度的一半
	{
		Complex Wn(cos(pi / mid), inv * sin(pi / mid));//单位根
		for (int r = mid << 1, j = 0; j < limit; j += r) //R是区间的长度，j表示前已经到哪个位置了
		{
			Complex w(1, 0);
			for (int k = 0; k < mid; ++k, w = w * Wn) //枚举左半部分
			{
                //complex x = a[j + k], y = w * a[j + mid + k];//蝴蝶效应   
                //a[j + k] = x+y; a[j + k + mid] = x-y; 
                //没有后效性也可以直接赋值给a，快一点
				buf[j + k] = a[j + k] + w * a[j + k + mid];
				buf[j + k + mid] = a[j + k] - w * a[j + k + mid]; 
			}
		}
		for (int i = 0; i < limit; ++i) a[i] = buf[i];
	}
}
signed main()
{
	n = read(), m = read();
	for (int i = 0; i <= n; ++i) a[i].x = read();
	for (int i = 0; i <= m; ++i) b[i].x = read();
	while (limit <= n + m) limit <<= 1, ++len;
	for (int i = 0; i < limit; ++i)
		pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (len - 1));
     // 在原序列中 i 与 i/2 的关系是 ： i可以看做是i/2的二进制上的每一位左移一位得来
    // 那么在反转后的数组中就需要右移一位，同时特殊处理一下奇数
	fft(a, 1);
	fft(b, 1);
	for (int i = 0; i <= limit; ++i) a[i] = a[i] * b[i];
	fft(a, -1);
	for (int i = 0; i <= n + m; ++i) 
        printf("%d ", (int)(a[i].x / limit + 0.5));
}

三次变两次优化

观察到上面的代码我们跑了三次肥肥兔$FFT$，现在我们有一种方法可以少跑一次

假设我们求$f(x)\ast g(x)$

设复多项式$h(x)=f(x)+i\ast g(x)$，$实部为f(x)，虚部为g(x)$

那么$h(x{)}^2=(f(x)+i\ast g(x){)}^2=f(x{)}^2-g(x{)}^2+i\ast 2\ast f(x)\ast g(x)$

我们只要把$h(x{)}^2$的虚部除以$2$就得到了结果

完全版：

注意三次变两次优化会令精度误差平方，请根据题目值域考虑是否使用

#include<bits/stdc++.h>
using namespace std;
#define eps (1e-8)
inline int read()
{
	int x = 0;
	char ch, f = 1;
	for (ch = getchar(); (ch < '0' || ch > '9') && ch != '-'; ch = getchar());
	if (ch == '-') f = 0, ch = getchar();
	while (ch >= '0' && ch <= '9')
	{
		x = (x << 1) + (x << 3) + ch - '0';
		ch = getchar();
	}
	return f ? x : -x;
}
const int N = 5e6 + 10;
const double pi = acos(-1.0);
int n, m;
int limit = 1, len;
int pos[N];
struct Complex
{
	double x, y;
	Complex(double tx = 0, double ty = 0)
	{
		x = tx, y = ty;
	}
	inline Complex operator + (const Complex t) const
	{
		return Complex(x + t.x, y + t.y);
	}
	inline Complex operator - (const Complex t) const
	{
		return Complex(x - t.x, y - t.y);
	}
	inline Complex operator * (const Complex t) const
	{
		return Complex(x * t.x - y * t.y, x * t.y + y * t.x);
	}
} a[N];
inline void fft(Complex *a, int inv)
{
	for (int i = 0; i < limit; ++i)
		if (i < pos[i]) swap(a[i], a[pos[i]]);//求出要迭代的序列
	for (int mid = 1; mid < limit; mid <<= 1)//待合并区间的长度的一半
	{
		Complex Wn(cos(pi / mid), inv * sin(pi / mid));//单位根
		for (int r = mid << 1, j = 0; j < limit; j += r)
		{
			Complex w(1, 0);
			for (int k = 0; k < mid; ++k, w = w * Wn)//枚举左半部分
			{
				Complex x = a[j + k], y = w * a[j + k + mid];//蝴蝶效应 
				a[j + k] = x + y;
				a[j + k + mid] = x - y;
			}
		}
	}
}
signed main()
{
	n = read(), m = read();
	for (int i = 0; i <= n; ++i) a[i].x = read();
	for (int i = 0; i <= m; ++i) a[i].y = read();
	while (limit <= n + m) limit <<= 1, ++len;
	for (int i = 0; i < limit; ++i)
		pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (len - 1));
	fft(a, 1);//三次变两次优化
	for (int i = 0; i <= limit; ++i) a[i] = a[i] * a[i];
	fft(a, -1);//三次变两次优化 
	for (int i = 0; i <= n + m; ++i) printf("%d ", (int)(a[i].y / limit / 2 + 0.5));
}

参考博客

%%%red

%%%自为风月马前卒