【FFT实战篇】C++实现：利用快速傅里叶变换快速计算（多项式）乘法

本文使用C++语言实现了快速傅里叶变换FFT并运用其进行多项式乘法计算，适合供学习过《数字信号处理》的读者学习使用FFT快速计算乘法。当然也可供对快速乘法感兴趣的没有系统学习过数字信号处理（从CFT到DTFT到DFT）的读者参考，提供一种不同于纯数学的推导思路。
另外本文着重探讨了实现二进制倒位序(bit-reversal)的方法。

快速傅里叶变换概念

快速傅里叶变换(FFT)是离散傅里叶变换(DFT)的快速算法实现，得以加速DFT是因为使用了分治的思想，本文的讲解采用典型的库利-图基(Cooley-Turkey)算法，即以时间抽选的FFT算法(DIT-FFT)。

但其实我们要关注的是FFT实现的目的：
加速离散傅里叶变换DFT($O(n^2)$)的计算，FFT的时间复杂度为$O(nlogn)$。

其中DFT中的离散是在离散时间傅里叶变换（DTFT）的基础上对频域（$2\pi$，周期）幅值进行离散得到的离散傅里叶级数DFS的基础上得到的（由于DFS时域频域都被一个周期内截断得到DFT，只考虑$0\sim N-1$的值，但也隐含周期性）所以，得到的是数字频率$\frac{2\pi }{N}$，$N$即为DFT的点数，在时域上也有重要性质，具体可查阅系列博客。

数字频率体现在表达式上就是$$X(k)=\sum _{n=0}^{N-1}x(n)W_N^{nk},k=0,1,2,\cdots ,N-1$$
这也就是DFT的表达式，其中$W_N^{nk}=e^{-j\frac{2\pi }{N}k}$是一个简略记号
可以看到离散的频率分量$k$也被限制在了一个周期范围内

同样的也有逆变换IDFT
$$x(n)=IDFT[X(k)]=\frac{1}{N}\sum _{n=0}^{N-1}X(k)W_N^{-nk},n=0,1,...,N-1$$

通过DFT我们就将离散序列从时域形式转化到了频域形式，这个频域由$2\pi$被$N$点抽样，目标就是求出周期内所有抽样点处的频域分量$X(k)$

想要了解原理，可参考【FFT理论篇】

快速傅里叶变换意义

通过$N$点DFT，将序列转化成了$N$个频域点表示。而通过IDFT，又可以将频域点重新转化成序列。
如果对卷积稍有了解的话应该知道连续情况下的卷积定理，时域上的卷积相当于频域上的相乘。
在DFT的情况下引入了圆周卷积和，情况稍有改变，但需要知道的是在一定条件下可以由圆周卷积和求出线性卷积和。

与多项式的关系

多项式$P(x)={\sum }_{n=0}^{N-1}{a}_i{x}^n$可由其系数向量$[a_{N-1},a_{N-2},\cdots ,a_0]$表示
两多项式相乘可通过系数向量相乘得出结果多项式的系数进而确定多项式

考虑$P(x)$与$Q(x)={\sum }_{n=0}^{M-1}{b}_i{x}^n$做多项式乘法
结果多项式的度数为$N+M-1$，对多项式系数进行乘法，相当于竖式乘法（每一项与另外一式中的所有系数相乘，最后相同幂次的结果相加），发现这其实就是计算卷积和中的“对位相乘相加法”[2017，程佩青，数字信号处理教程，清华大学出版社]，因此竖式乘法/多项式乘法可以当做卷积和运算来处理。

下面是一个例子，通过图1（下方)更好地去体会

结论一：多项式乘法并不是按照幂的次数来两两相乘多项式系数的数乘，而是多项式系数序列的卷积和

在时域上计算卷积和需要$O(n^2)$的时间复杂度，因此考虑利用卷积定理将其转化到频域($O(nlogn)$)上的数乘进行计算。也就是先对两个多项式系数序列进行FFT转成频域表示，然后将FFT的结果进行$O(n)$的相乘，再将数乘结果作IFFT转成时域表示得到了多项式乘法的结果。

离散傅里叶变换的圆周卷积和定理

既然目标是计算两个多项式系数序列的卷积（线性卷积和），现在关注如何通过FFT来快速计算。
引出DFT的圆周卷积和定理

长度为$N_1$的序列$x_1(n)$与长度为$N_2$的序列$x_2(n)$的$L$点圆周卷积和为
$\begin{array}{rl}y(n)=& x_1(n)Ⓛx_2(n)=x_2(n)Ⓛx_1(n)\\ =& [\sum _{m=0}^{l-1}{x}_1(m)x_2((n-m){)}_L]R_L(n),L\ge max[N_1,N_2]\end{array}$

$L$为圆周卷积和的点数。通过表达式看出圆周卷积和计算公式与DFT的形式类似，都是在有限点数内进行截断得到的结果。

之所以与连续卷积情况不同是因为DFT本身是在一个周期内截断并且具有周期性的结果，产生了与圆周（周期）有关的系列性质。

我们要求的线性卷积和的长度为$N_1+N_2-1$，只有当取$L\ge N_1+N_2-1$时，才能从圆周中恢复出完整的线性卷积和序列。为了满足FFT要求的点数，我们要使DFT（也就是FFT）的点数，同时也是圆周卷积和的点数$N=L=2^m\ge N_1+N_2-1,m为正整数$

圆周卷积和定理

$$y(n)=x_1(n)Ⓛx_2(n)$$
$$Y(k)=X_1(k)X_2(k),L点$$
上面我们已经知道了由圆周卷积和可以得到我们要求的线性卷积和，根据圆周卷积定理公式我们得到第二个结论

结论二：为了计算线性卷积和，需要将两个序列补到一个为2的整次幂的点数$N=2^m\ge N_1+N_2-1,$并对两个序列作$N$点DFT，以通过圆周卷积和定理计算$N$点圆周卷积。对$N$点圆周卷积作$N_1+N_2-1$点截断就是线性卷积和。

情境分析

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

题目描述

给定一个 $n$ 次多项式 $F(x)$，和一个 $m$ 次多项式 $G(x)$。

请求出 $F(x)$ 和 $G(x)$ 的卷积。

输入格式

第一行两个整数 $n,m$。

接下来一行 $n+1$ 个数字，从低到高表示 $F(x)$ 的系数。

接下来一行 $m+1$ 个数字，从低到高表示 $G(x)$ 的系数。

输出格式

一行 $n+m+1$ 个数字，从低到高表示 $F(x)\cdot G(x)$ 的系数。

样例 #1

样例输入 #1

1 2
1 2
1 2 1

样例输出 #1

1 4 5 2

提示

保证输入中的系数大于等于 $0$ 且小于等于 $9$。

对于 $100\%$ 的数据：$1\le n,m\le {10}^6$。

分析

按照上面FFT与多项式原理的介绍，将输入的两个多项式序列作FFT，再对FFT的结果数乘后作IFFT得到结果。
$$\begin{array}{rl}X(k)=& DFT[x(n)]=\begin{array}{c}\underbrace{\sum _{n=0}^{N-1}x(n)W_N^{nk}}\\ N点长序列\end{array}=\begin{array}{c}\underbrace{\sum _{n=0}^{N-1}x(n)W_N^{nk}}\\ n为偶数\end{array}+\begin{array}{c}\underbrace{\sum _{n=0}^{N-1}x(n)W_N^{nk}}\\ n为奇数\end{array}\\ =& \sum _{r=0}^{\frac{N}{2}-1}x(2r)W_N^{2rk}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_N^{2rk}\\ =& \sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_N^{r{k}_1}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_N^{r{k}_1},k_1=2k=0,2,4,\cdots ,2(N-1)\\ =& \begin{array}{c}\underbrace{\sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{rk}}\\ \frac{N}{2}点长序列\end{array}+W_N^k\begin{array}{c}\underbrace{\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{rk}}\\ \frac{N}{2}点长序列\end{array},k=0,1,2,\cdots ,N-1\\ =& \{\begin{array}{l}{\sum }_{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{rk}+W_N^k{\sum }_{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{rk},k=0,1,2,\cdots ,\frac{N}{2}-1\\ {\sum }_{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{rk}-W_N^k{\sum }_{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{rk},k=\frac{N}{2},\frac{N}{2}+1,\cdots ,N-1\end{array}\end{array}$$

Radix-2 DIT-FFT的思路就是将计算$N$点DFT分解成计算两个$\frac{N}{2}$点DFT之和，如此分解进行下去直到…计算出2点DFT（不可再分解)，再按照最后的公式一层层地往回代出结果。
这样子一层层将一次$O(n^2)$的运算分解成缩小的$n$的两次运算，将大问题分解成了子问题，由于进行到底，每一层只有$O(n)$的运算，这也是分治法(Divide And Conquer)的基本思路。

分解的原理是按照点数（时间）的奇偶对DFT计算式进行重排，$N$点序列的$N$点DFT变成了两个$\frac{N}{2}$点长的序列的$N$点DFT，由于$\frac{N}{2}$点DFT只有$\frac{N}{2}$个不重复的频率，推导的最后变成了计算两个$\frac{N}{2}$点序列的$\frac{N}{2}$点DFT，但是要分段计算。
要分段是因为按照时间抽选，从奇数序列DFT的$x(2r+1)W_N^{(2r+1)k}$中提出$W_N^k$，留下${\sum }_{r=0}^{\frac{N}{2}-1}x(2r+1)W_N^{2rk}$使得频率分量与偶数序列保持一致。
这个提出来的复数系数$W_N^k=e^{-j\frac{2\pi }{N}k}$有时也被称为旋转因子(twiddle factor)
而分段是用到了$W_N^k$的一个性质：$W_{N/2}^{k+\frac{N}{2}}=-W_{N/2}^k$，$\frac{N}{2}$之后要变号

在上面完整的推导过程中还用到了$W_N^{nk}$的可约性：$W_N^{in}=W_{\frac{N}{i}}^n$，$W_{iN}^{in}=W_N^n$

这里顺便注意一个点，如果在第一段计算旋转因子时记录下$W_N^k{\sum }_{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{rk}$的值，在计算第二段时可以直接使用，省去了一次乘法。${\sum }_{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{rk}$同理。

这样下来分解一次DFT只需要进行两次复数加法，一次复数乘法。（注意${\sum }_{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{rk}$和${\sum }_{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{rk}$是DFT的结果由下一层在返回时已经给出为已知量，实际上是一次复数加法，一次复数减法，再加上旋转因子乘上$x_2(r)$DFT结果的复数乘法）。
这一过程若画出蝴蝶图(butterfly diagram)能更好的看出特点，对于整体复杂度的分析也更容易把握，蝴蝶图请参阅理论篇。其实蝴蝶图就是横版的图2，理解蝴蝶运算需要把握三点：
一，序列的总长度不变，每一层的序列总点数是$N$，且蝴蝶运算结的数量每层相等
二，序列的层数为$lo{g}_{2}N$，这也是对应的下标二进制表示的位数
三，计算每一层的“父节点”的值只在上一层的下标的“中点“($\frac{N}{2}$)处考虑分段，且当前层的两个节点各长$\frac{N}{2}$要各自对父节点做出贡献，在图中这两个节点的元素之间的距离是随层数规律变化的

由于涉及到复频率，代码中需要使用复数，这里使用C++ STL中的complex，具体使用方法见代码

递归版FFT

在递归函数中直接开辟出新的两个数组保存重排结果并作为参数传递下去
旋转因子可以在循环中不断累乘得到

void FFT(complex<double> A[], int n) 
{
    if (n == 1)
        return;
    complex<double> AE[n / 2], AO[n / 2];
    for (int i = 0; i < n / 2; i++)
    {
        AE[i] = A[i * 2]; //偶数频率，这里奇偶指的是对于当前序列
        AO[i] = A[i * 2 + 1]; //奇数频率
    }
    FFT(AE, n / 2);
    FFT(AO, n / 2); //递归求解
    complex<double> W, Wk;
    W.real(cos(pi * 2 / n));
    W.imag(-sin(pi * 2 / n));
    Wk = 1; //旋转因子初始值
    for (int i = 0; i < n / 2; i++)
    {
        A[i] = AE[i] + Wk * AO[i];
        A[i + n / 2] = AE[i] - Wk * AO[i];
        Wk *= W;
    }
}

上面的代码比较朴素，实际上可以进行很多优化。
最致命的是这段代码的内存开销过大，在本题中$N=1e6$下在递归函数下（栈区）开辟一个$5e5$大小的complex甚至会直接爆栈，为此我们引入迭代版FFT。

但是想要迭代循环实现FFT，我们有一个至关重要的问题函待解决。我们之所以进行递归，就是为了不断分解序列（二叉树的后序遍历）得到最底下一层的序列的排序（它们对应在根节点的排列下的下标），我们把这个顺序叫做倒位序，因为通过观察发现其二进制表示为自然序二进制的倒位。

为什么是自然序二进制的倒位呢，我们在每一层向下分解时实际上是按照奇偶进行分组，也就是以二进制的0和1作为父节点不断向下构建二叉树，从图4的下标生成树可以清晰理解二进制倒位序背后的原因

其实在自顶向下得到这个倒位序的过程中并没有进行运算，我们只需关心从倒位序向上合并得到结果的过程，所以要找到一种能够直接求出倒位序的方法。

二进制倒位序的算法

我们需要在原址（即输入序列的数组上）计算出其倒位序。考虑朴素的算法，使用递归或者循环，依次从右向左访问二进制的每一位并将其保存或与对称的对应位交换，这样要计算出从0到N-1的每个数的倒位序的时间复杂度是$O(N^2)$，速度太慢。
这里介绍几种容易理解但又比朴素方法高效的二进制倒位序算法。

Gold-Rader算法

雷德算法(Gold-Rader‘s algorithm)是最早提出也是运用最广泛的计算倒位序的方法[1969,B,Gold and C.M.Rader,Digital Processing of Signals,McGraw-Hill]，广大《数字信号处理》教材对这种方法有初步的介绍。
雷德算法相较于朴素算法的改进为，不再依次正序枚举所有要进行倒位的数再逐位反转，而是考虑观察规律直接按序枚举倒位序数。观察倒位序数的排列发现，倒位序的产生规律与自然序数类似，只是进位方向相反，雷德算法模拟倒位序的进位规则产生下一个倒位序数。

由图5观察到一些进位规律（在十进制下）：最高位对应$\frac{N}{2}$，若最高位为0，优先进最高位；若最高位为1，最高位归零去进下一位（直到从右向左找到不为一的首个零位，遇到的1全部归零），次高位对应$\frac{N}{4}$，然后是$\frac{N}{8}$，…，依次类推
下面给出雷德算法的实现

int Rader(int num,int len) //Rader算法产生num的下一个倒位序数
{
    int k = len >> 1; //len为DFT的点数，是十进制下的最大值加一
    while(num >= k){ //第一种情况，进位后当前位需要置零
        num -= k;![请添加图片描述](https://img-blog.csdnimg.cn/direct/18cecc24070a4231a9aca726d477ddfa.png)

        k >>= 1;
    }
    if(num < k)num += k; //第二种情况，只需要最高位进位
    return num;
}

倒位序序列与二进制位数的关系

雷德算法找到了相邻两个倒位序二进制数的联系，通过模拟进位的方法可以顺序求出下一个数，但是在一些情况下仍然需要遍历几乎全部数位。
其实不同长度之间二进制倒序序列是有联系的，现在考虑找出规律，由当前长度的倒序序列推出当前二进制位长度加一的倒序序列。

规律及证明

考虑二进制位的长度为$M$，则十进制序列的长度（点数）$N=2^M$，列出倒序序列随$M$变化的表格如下

其实这就是图4下标生成树的表格形式，我们容易看出：
十进制点长$N=2^M$ 的倒序序列(以下称为目标序列)可由点长为$N_1=\frac{N}{2}$（二进制长度为$M-1$）的倒序序列(以下称为生成序列)按照如下规则生成（十进制视角）：
1.生成序列的每个元素$n$ 不变
2.同时在每个元素$n$后面插入一个新的元素$n+N_1$

刘大庆等人证明了这条性质[2018.刘大庆，林浩然，陈树越.快速傅里叶变换中计算倒序的新思路[J].电子与信息学报]

证明的简要思路如下：
1.将目标序列（二进制长度为$M$)的$N$个元素每相邻两个分成一组，共分成$N_1=\frac{N}{2}$组
2.第$n$组的元素为$a_{2n}$,$a_{2n+1}$，$n=0,1,2,\cdots ,\frac{N}{2}-1$，这里的下标是从左往右数的当前序列的下标
3.考虑十进制元素$n$，在目标序列长度下的二进制自然序表达为$(0,n_{M-2},n_{M-3},\cdots ,n_1,n_0)$（$n<\frac{N}{2}=2^M$)
4.则生成序列长度下中$n$可自然序表达为$(n_{M-2},n_{M-3},n_{M-4},\cdots ,n_1,n_0)$，倒序表达为$(n_0,n_1,n_2,\cdots ,n_{M-3},n_{M-2})$
5.倒位序序列中的第$n$个元素$a_n$根据定义恰为$n$的二进制倒位序，因为倒位序序列是自然序序列的倒序而自然序序列从$0$到$N-1$，$n$包含在内
6.根据第5点观察，生成序列的第$n$个元素$a_n$即为$(n_0,n_1,n_2,\cdots ,n_{M-3},n_{M-2})$
7.根据二进制移位的性质（移位后倒序），由生成序列的$a_n$移位得到的二进制$a_{2n}$即为$(0,n_0,n_1,n_2,\cdots ,n_{M-3},n_{M-2})$，这与$(n_0,n_1,n_2,\cdots ,n_{M-3},n_{M-2})$即目标序列下$a_{2n}$的值是一致的，证明了生成规则的第1点
8.先移位后按照规律进位，$a_{2n+1}$为$(1,n_0,n_1,n_2,\cdots ,n_{M-3},n_{M-2})$，这是生成序列中没有的，转化成十进制表示即证明了生成规则的第2点

Rius等人也发现了类似的规律，他们是从同一二进制长度下分组的思路来实现的。[1995.Juan M.Rius and R.De Porrata-Dòria.New FFT Bit-Reversal Algorithm[J].IEEE Transactions on Signal Processing]

再次改进及代码实现

上述的方法虽然发现了倒位序序列与二进制长度之间的联系，但是没有精确到元素，需要对整段序列中的所有元素进行操作，产生了大量不同长度的中间元素，占用额外的内存，第一篇文章中的流程图说明了这一点。
我们深入挖掘第一种方法背后的原理发现其实可以使用递推的方法依次求出倒位序数，整个过程不需要额外的空间，是真正$O(n)$的实现。

1.生成序列的每个元素$n$ 不变
2.同时在每个元素$n$后面插入一个新的元素$n+N_1$

目标序列$b$的十进制长度是生成序列$a$的两倍，由证明过程可知 $b_{2n}=a_n$,$b_{2n+1}=a_n+N_1$
十进制数$n$的一半$\frac{n}{2}$的二进制长度减一，因此这两个序列可以合二为一为$R$

在求倒位序第$n$个元素$R[n]$递推过程中已知$R[n/2]$即倒位序的第$\frac{n}{2}$个元素。若是第1种情况，直接维持元素不变；若为第二种情况，需要加上$N/2$

$$R[n]=R[n>>1]>>1+(R\&1)\ast n/2$$

这是一种非常巧妙的方法，只需根据倒位序的含义就可以理解，但是要详细掌握第一种方法的推导过程，详情请见下面的代码。

迭代版FFT

在原址上倒位序后仅需自底向上进行蝴蝶运算，注意原址运算的先后覆盖问题，代码中使用了额外的变量来处理这一问题。
在蝴蝶运算中使用了两个不同的指针，这部分编写有一定的技巧。

void FFT_Iteration(complex<double> A[], int n)
{
    // bit reversal
    pos[0] = 0;//第0个元素倒位序已知
    for (int i = 1; i < n; i++)
        pos[i] = pos[i / 2] / 2 + (i % 2) * n / 2; //求倒位序
    // cout << endl;
    for (int i = 0; i < n; i++)
        if (i < pos[i])
            swap(A[i], A[pos[i]]); //根据倒位序原址交换
    // butterfly diagram
    for (int l = 2; l <= n; l *= 2)
    {
        complex<double> W, Wk;
        for (int i = 0; i < n; i += l)
        {
            Wk = 1;
            W.real(cos(pi * 2 / l));
            W.imag(-sin(pi * 2 / l));
            for (int j = 0; j < l / 2; j++)
            {
                // A[i+j]=A[i+j]+Wk*A[i+j+l/2]; 错误写法
                // A[i+j+l/2]=A[i+j]-Wk*A[i+j+l/2];
                complex<double> x = A[i + j], y = A[i + j + l / 2] * Wk;//使用额外变量
                A[i + j] = x + y;
                A[i + j + l / 2] = x - y;
                Wk *= W;
            }
        }
    }
}

多项式乘法完整代码（迭代版）

IFFT可以直接按照公式，将FFT中的$W_N^k$换成$W_N^{-k}$，再在结果处乘上$\frac{1}{N}$即可。
或者为了保持IFFT函数的对应性，注意到$N=2^L$（$L$为FFT层数），在每层计算时结果累乘$\frac{1}{2}$也不失为一种好方法。
实际上计算两个实数序列（本题）的FFT还可以进行优化到只需一次FFT，这里篇幅受限不再介绍。

#include <bits/stdc++.h>
using namespace std;
const int N = 4e6 + 10;
const double pi = acos(-1);
#define endl '\n'
complex<double> A[N], B[N], C[N];
int pos[N];
void FFT_Iteration(complex<double> A[], int n)
{
    // bit reversal
    pos[0] = 0;
    // cout << "after reversal ";
    for (int i = 1; i < n; i++)
    {
        pos[i] = pos[i / 2] / 2 + (i % 2) * n / 2;
        // cout << pos[i] << " ";
    }
    // cout << endl;
    for (int i = 0; i < n; i++)
        if (i < pos[i])
            swap(A[i], A[pos[i]]);
    // for (int i = 0; i < n; i++)
    // cout << A[i] << " ";
    // cout << endl;
    // butterfly diagram
    for (int l = 2; l <= n; l *= 2)
    {
        complex<double> W, Wk;
        for (int i = 0; i < n; i += l)
        {
            Wk = 1;
            W.real(cos(pi * 2 / l));
            W.imag(-sin(pi * 2 / l));
            for (int j = 0; j < l / 2; j++)
            {
                // A[i+j]=A[i+j]+Wk*A[i+j+l/2];
                // A[i+j+l/2]=A[i+j]-Wk*A[i+j+l/2];
                complex<double> x = A[i + j], y = A[i + j + l / 2] * Wk;
                A[i + j] = x + y;
                A[i + j + l / 2] = x - y;
                Wk *= W;
            }
        }
    }
}
void IFFT_Iteration(complex<double> A[], int n)
{
    // bit reversal
    pos[0] = 0;
    // cout << "after reversal ";
    for (int i = 1; i < n; i++)
    {
        pos[i] = pos[i / 2] / 2 + (i % 2) * n / 2;
        // cout << pos[i] << " ";
    }
    // cout << endl;
    for (int i = 0; i < n; i++)
        if (i < pos[i])
            swap(A[i], A[pos[i]]);
    // for (int i = 0; i < n; i++)
    // cout << A[i] << " ";
    // cout << endl;
    // butterfly diagram
    for (int l = 2; l <= n; l *= 2)
    {
        complex<double> W, Wk;
        for (int i = 0; i < n; i += l)
        {
            Wk = 1;
            W.real(cos(pi * 2 / l));
            W.imag(sin(pi * 2 / l));
            for (int j = 0; j < l / 2; j++)
            {
                // A[i+j]=A[i+j]+Wk*A[i+j+l/2];
                // A[i+j+l/2]=A[i+j]-Wk*A[i+j+l/2];
                complex<double> x = A[i + j], y = A[i + j + l / 2] * Wk;
                A[i + j] = x + y;
                A[i + j + l / 2] = x - y;
                Wk *= W;
            }
        }
    }
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    int n, m;
    cin >> n >> m;
    for (int i = 0; i <= n; i++)
        cin >> A[i];
    for (int i = 0; i <= m; i++)
        cin >> B[i];
    int N = 1;
    while (N <= n + m + 1)
        N *= 2;
    FFT_Iteration(A, N);
    FFT_Iteration(B, N);
    for (int i = 0; i < N; i++)
        C[i] = A[i] * B[i];
    IFFT_Iteration(C, N);
    for (int i = 0; i < n + m + 1; i++)
        cout << int(C[i].real() / N + 0.5) << " ";
    // cout << round(C[i].real() / N ) << " "; 速度慢，且有-0问题
    return 0;
}