[FFT/IFFT]快速傅里叶(逆)变化 + 递归和递推模板

现在时间是2021-2-2，重新回来看2019学习的一知半解的$FFT$，又有了新的理解
所以修改了以往写过的文章，并增添些许内容
因为过去一年多，上了高中，学的知识多了些，以前不懂的有些东西现在看来挺简单的？？

• Add
  建议了解系数和点值表示法后直接从复数看起
  因为前面很多是第一次学习的时候，为了全面了解
  然而似乎并没有起到这个效果？？

引入

如果给你一个多项式 $A(x)=\sum a_n{x}^n$ 和 $B(x)=\sum b_n{x}^n$，求 $A(x)\cdot B(x)$
你会怎么做？？
——可能只能选择$O(n^2)$，做${\sum }_{i=1}^n{\sum }_{j=1}^{m}ai\ast bj$
但是你觉得那些毒瘤出题dalao，会让你轻轻松松$O(n^2)$水过去吗？

如此之高的时间复杂度永远成为了多项式乘法的一个瓶颈…
直到伟大的$FFT$就出现了，将其优化到$O(nlogn)$

---

系数和点值表示法

对于求一个$n-1$次的多项式$f(x)$，可以有两种表示方法，并且可以互相转化

系数表示法

概念

$$f(x)=\sum _{i=0}^{n-1}{c}_i\ast x^i$$
什么意思呢？
🍔设$f(x)=(a_0\ast x^0+a_1\ast x^1+a_2\ast x^2)\ast (b_0\ast x^0+b_1\ast x^1+b_2\ast x^2)$

那么一一相乘将之拆成了九项式子相加
$f(x)=a_0{b}_0\ast x^0+a_0{b}_1\ast x^1+a_0{b}_2\ast x^2+a_1{b}_0\ast x^1+a_1{b}_1\ast x^2$
$+a_1{b}_2\ast x^3+a_2{b}_0\ast x^2+a_2{b}_1\ast x^3+a_2{b}_2\ast x^4$

进行同类项合并
$f(x)=a_0{b}_0\ast x^0+(a_0{b}_1+a_1{b}_0)\ast x_1+(a_0{b}_2+a_1{b}_1+a_2{b}_0)\ast x^2$
$+(a_1{b}_2+a_2{b}_1)\ast x^3+a_2{b}_2\ast x^4$

$c_0=a_0{b}_0$
$c_1=(a_0{b}_1+a_1{b}_0)$
$c_2=(a_0{b}_2+a_1{b}_1+a_2{b}_0)$
$c_3=(a_1{b}_2+a_2{b}_1)$
$c_4=a_2{b}_2$

用中国话来理解就是$c_i$是我们整合后的对于$x^i$的系数和，将这些相加就是最后的$f(x)$

---

• Add
  $c_i$本质来说可以看作是两个函数的卷积
  $$h(n)=\sum _{i=0}^{n-1}{c}_i{x}^i,f(n)=\sum _{i=0}^{n-1}{a}_i{x}^i,g(n)=\sum _{i=0}^{n-1}{b}_i{x}^i$$$$\Rightarrow c_i=\sum _{j=0}^i{a}_j\times b_{i-j}$$

---

系数表示法$\to$点值表示法

就是把$x_i$带进去就可以算出来每一个点$i$的函数值，就可以表示该点$(x_i,y_i)$
$$y_i=\sum _{j=0}^{n-1}{c}_j\ast x_i^j$$

系数表示法的优缺点

优点
多项式的求值计算效率高，对于$A(x)=a_0+a_1\ast x^1+a_2\ast x^2...a_n\ast x^n$
提一个同类项$x$，变成
$A(x)=a_0+x(a_1+a_2\ast x+...a_n\ast x^{n-1})$，不停地提$x$出来
我们就可以在$a_0$处通过霍纳法则$O(n)$算出来

多项式的加减计算效率也高
$A(x)=a_0+a_1\ast x^1+a_2\ast x^2+...+a_{n-1}\ast x^{n-1}$
$B(x)=b_0+b_1\ast x^1+...+b_{n-1}\ast x^{n-1}$
可以通过$O(n)$，算出$C(x)=c_0+c_1\ast x^1+...+c_{n-1}\ast x^{n-1}$
对于每一个$i\in [0,n)$，都有$c_i=a_i+b_i$，其实就是直接系数方面的相加减

缺点
多项式的乘法计算时间复杂度将达到$O(n^2)$
1.感性理解就是我们要枚举$A$里面的每一项，再与$B$里面的每一项进行相乘再合并同类项
2.数学公式表达则是：解释一下为什么上界是$2n-2$
额其实很好想，$A,B$的最高项都是$x^{n-1}$相乘肯定就是$C$的最高项，也就是$x^{2n-2}$
$$C(x)=\sum _{i=0}^{2n-2}{c}_i\ast x^i,c_i=\sum _{j=0}^i{a}_j{b}_{i-j}$$

点值表示法

概念

给一堆点对$(x_1,x_2,x_3...x_n),(y_1,y_2,y_3...y_n)$，满足$f(x_i)=y_i$
即$(x_i,y_i)$是曲线上$y=f(x)$的点

---

• Add
  这样表示似乎更好？
  $$(x_0,y_0),(x_1,y_1),(x_2,y_2)....(x_{n-1},y_{n-1})$$

---

用中国话来讲就是我们知道了平面直角坐标系上某条函数的$n$对点
然后就可以勾勒出这一条唯一的函数图象

要确定一条函数的图像，要至少知道函数最高次+1个不同的点

简单证明一下：
1.感性理解，我们说两个点确定一条直线，也就是说要两个点才能画出一次函数
而我们的抛物线又要三个点才能画出二次函数$...$以此类推
2.运用解方程的方法，我们面对四个点会设$f(x)=a\ast x^3+b\ast x^2+c\ast x^1+d$
四个不同的方程对应四个解
感觉好像是一样的证明，别管这么多了，反正都是简单证明，口胡口胡

点值表达式–>系数表达式

$$f(x)=\sum _{i=0}^{n-1}{y}_i\frac{{\prod }_{j\ne i}(x-x_j)}{{\prod }_{j\ne i}(x_i-x_j)}$$
这个证明要用到拉格朗日插值法，但是因为我们一般不用这玩意儿，老子就不搞了，太现实了

点值表达式的优缺点

优点
加减法计算效率高：对两个点值表达的次数界为$n$的多项式，计算只有$O(n)$
如果$C(x)=A(x)+b(x)$，那么$C(x_k)=A(x_k)+B(x_k)$

更具体而言：对于给定的$A(x_0,y_0),(x_1,y_1)...(x_{n-1},y_{n-1})),B(x_0,y_0^{\prime }),(x_1,y_1^{\prime })...(x_{n-1},y_{n-1}^{\prime }))$
那么$A$和$B$对相同的$n$个点对求和，$C$的点对就表示成$C(x_0,y_0+y_0^{\prime }),(x_1,y_1+y_1^{\prime })...(x_{n-1},y_{n-1}+y_{n-1}^{\prime }))$

乘法计算效率也高：对两个点值表达的次数界为$n$的多项式，计算只有$O(n)$
如果$C(x)=A(x)\ast b(x)$，那么$C(x_k)=A(x_k)\ast B(x_k)$
这样只需要将$A,B$进行逐点相乘就可以求出了$C$，但是$C$的次数界要达到$2n$
$A,B$次数界也只有$n$，所以我们必须对$A,B$进行扩点处理，将其扩大成$C$的次数界

更具体而言：扩充$A(x_0,y_0),(x_1,y_1)...(x_{2n-1},y_{2n-1})),B(x_0,y_0^{\prime }),(x_1,y_1^{\prime })...(x_{2n-1},y_{2n-1}^{\prime }))$
$C$的点对就表示成$C(x_0,y_0+y_0^{\prime }),(x_1,y_1+y_1^{\prime })...(x_{n-1},y_{2n-1}+y_{2n-1}^{\prime }))$

缺点
我们如何求一个新点的值呢？是不是只能转化成系数表达式，用$O(n)$计算
但是时间复杂度就在转化这里达到了$O(n^2)$

换言之：对于多项式 $A(x)$ 和 $B(x)$，假设 $degA+degB<n$（
deg是数学中的表示多项式的次数的玩意儿）
如果有 $A$ 和 $B$ 在 $x_0,x_1,...,x_{n-1}$ 处的点值表示
则 $(A\cdot B)$的点值表示可以通过$(A\cdot B)(x_i)=A(x_i)\cdot B(x_i)$ 在 $O(N)$ 时间内得到
还原 $(A\cdot B)$ 为系数表示就实现了多项式乘法，但是还原的时间$O(n^2)$

🍔：所以如果有一道题给我们系数表达式，最后又让我们输出结果的系数表达式
我们用以上的方法虽然计算成点值表达式只用了$O(n)$
但是最后在转化成系数表达式的时候，时间复杂度还是蹭蹭蹭地涨到了$O(n^2)$
看上面的式子，我们会面临枚举$i,j$的难题，还是没有在本质上解决问题

---

但是我们的$FFT$就剋以做到以上的转化且只用$O(nlogn)$

1.把已知的一个多项式转化成对应的点值表示
2.把已知的点值表示转换成对应的多项式

说了这么多，还是没有告诉我怎么做啊！！！不急慢慢往下看

---

复数和单位复根

复数

我们把形如$z=a+bi$（a,b均为实数）的数称为复数，其中$a$称为实部，$b$称为虚部，$i$称为虚数单位
$i=\sqrt{-1}$，即$i^2=-1$

我们可以把复数当做一个向量丢在二维平面，即平面直角坐标系

---

百度百科说：复数之间的加减乘（除）是可以直接算的，除法因为不怎么用就不说了

1.加法法则
设$z1=a+bi，z2=c+di$是任意两个复数，
则它们的和是 $(a+bi)+(c+di)=(a+c)+(b+d)i$
两个复数的和依然是复数，它的实部是原来两个复数实部的和，它的虚部是原来两个虚部的和，当然复数的加法满足交换律和结合律

---

2.减法法则
设$z1=a+bi，z2=c+di$是任意两个复数，
则它们的差是 $(a+bi)-(c+di)=(a-c)+(b-d)i$
两个复数的差依然是复数，它的实部是原来两个复数实部的差，它的虚部是原来两个虚部的差

---

3.乘法法则
设$z1=a+bi，z2=c+di(a、b、c、d\in R)$是任意两个复数，
那么它们的积$(a+bi)(c+di)=(ac-bd)+(bc+ad)i$
其实就是把两个复数相乘，类似两个多项式相乘，展开得: $ac+adi+bci+bd{i}^2$，因为$i^2=-1$，所以结果是$(ac－bd)+(bc+ad)i$ ，两个复数的积仍然是一个复数
此时，复数相乘表现为幅角相加，模长相乘

单位复根

定义$n$次单位复根为${\omega }_n^i$，满足$x^n=1$的复数$x$，表现在平面直角坐标系中

单位复根满足的性质如下：可以想象成在单位圆上的旋转

${\omega }_n^a\ast {\omega }_n^b={\omega }_n^{a+b}$
${\omega }_n^i={\omega }_n^{i+n}$，也就是转了一圈单位圆最后在坐标系上只转了幅角为i
${\omega }_{kn}^{ki}={\omega }_n^i$
${\omega }_n^i=-{\omega }_n^{i+n/2}$，$-$可以理解为倒着转

因为单位复根刚好有$n$个，可以分别一一对应我们的$n-1$次多项式，形成点值表达式

{ ( ω n 0 , f ( ω n 0 ) ) , ( ω n 1 , f ( ω n 1 ) ) , ( ω n 2 , f ( ω n 2 ) ) . . . , ( ω n n − 1 , f ( ω n n − 1 ) ) } \{(\omega_n^0,f(\omega_n^0)),(\omega_n^1,f(\omega_n^1)),(\omega_n^2,f(\omega_n^2))...,(\omega_n^{n-1},f(\omega_n^{n-1}))\} {(ωn0​,f(ωn0​)),(ωn1​,f(ωn1​)),(ωn2​,f(ωn2​))...,(ωnn−1​,f(ωnn−1​))}

我们的$FFT$是要求$n$为2的幂次的

---

• Add
  对于一个函数$f(n)={\sum }_{i=0}^{n-1}{c}_i{x}_i$
  其实可以上调$n$但不能下降$n$
  意思就是可以将$n$配成任何比$n$大的$n^{\prime }$，系数$c$直接配成$0$，不就行了？
  所以$FFT$的$n$的要求是完全可以人为调控达到的

---

所以像上图的五个单位复根
其实我们是分成了八个单位复根，然后就只用前五个
如图分成了八份，只用其中涂绿了的五份

---

• Add
  前面提到是将复数当成向量放在二维平面的单位圆里
  所以对于单位圆上的一点，假设角度为$x$，那么该点可以表示为$(cosx,isinx)$
  对于两个角度分别为$x,y$的在单位圆上的点，相乘
  $(cosx,isinx)\times (cosy,isiny)=$
  $(cosxcosy-sinxsiny,i(sinxcosy+cosxsiny))$
  发现就是角度为$x+y$的点坐标$(cos(x+y),isin(x+y))$
  这也恰恰应证了复数相乘表现为幅角相加，模长相乘
  点乘算出来的结果是一个点
  叉乘算出来的结果仍是一个向量

---

傅里叶正变换（一般形式$\to$点值表达式）

理论

---

• Add
  $FFT$就是知道用点值表达式表示函数

---

FFT 的正变换实现，是基于对多项式进行奇偶项分开递归再合并的分治进行的
对于 $n-1$ 次多项式，我们选择插入 $n$ 次单位根求出其点值表达式，
这就跟我们引入单位复根的原因相结合了

设$f(x)=a_0+a_1\ast x^1+a_2\ast x^2+...+a_{n-1}\ast x^{n-1}$，
$f_0(x)=a_0+a_2\ast x+a_4\ast x^2+a_6\ast x^3+...$，$f_1(x)=a_1+a_3\ast x+a_5\ast x^2+...$
则$f(x)=f_0(x^2)+x\ast f_1(x^2)$
证明的话把这个式子展开就行了，跳过

也就是说我们把$f(x)$分成了两类，奇数项分成一类，偶数项分成一类，去得到上列公式

接着，令$n=2\ast p$，那么就有以下转化
$$f({\omega }_n^i)=f_0(({\omega }_{n/2}^{i/2}{)}^2)+{\omega }_n^i\ast f1((w_{n/2}^{i/2}{)}^2)=f_0({\omega }_p^i)+{\omega }_n^i\ast f_1({\omega }_p^i)$$
$$f({\omega }_n^{i+p})=f_0(({\omega }_{n/2}^{(i+p)/2}{)}^2)+{\omega }_n^{i+p}\ast f1((w_{n/2}^{(i+p)/2}{)}^2)$$$$=f_0({\omega }_p^{i+p})+{\omega }_n^{i+p}\ast f_1({\omega }_p^{i+p})=f_0({\omega }_p^i)-{\omega }_n^i\ast f_1({\omega }_p^i)$$

由上述式子我们可以知道，如果我们知道$f_0({\omega }_p^i),f_1({\omega }_p^i)$，我们就可以$O(1)$算出$f({\omega }_n^i),f({\omega }_n^{i+n/2})$
那么如果我们递归求出了$f_0(x),f_1(x)$的$n/2$次的插值，我们就能$O(n)$的算出$f(x)$的$n$次单位根的插值，所以时间复杂度则是$O(nlogn)$

一言以蔽之：当 $x$ 取遍所有 $n$ 次单位复根时，$x^2$ 取遍所有 $(n/2)$ 次单位复根

递归模板

struct complex {//先自己手打STL里面的复数，可以防止某些**卡常
	double real, i;
	complex () {}
	complex ( double xx, double yy ) {
		real = xx;
		i = yy;
	}
}a[MAXN], b[MAXN];

complex operator + ( complex s, complex t ) {
	return complex ( s.real + t.real, s.i + t.i );
}
complex operator - ( complex s, complex t ) {
	return complex ( s.real - t.real, s.i - t.i );
}
complex operator * ( complex s, complex t ) {
	return complex ( s.real * t.real - s.i * t.i, s.real * t.i + s.i * t.real );
}

const double pi = acos ( -1.0 );

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 w = complex ( cos ( 2 * pi / limit ), inv * sin ( 2 * pi / limit ) ), p = complex ( 1, 0 );
	for ( int i = 0;i < ( limit >> 1 );i ++, p = p * w ) {
		a[i] = a1[i] + p * a2[i];
		a[i + ( limit >> 1)] = a1[i] - p * a2[i];
	}
}

---

傅里叶逆变换（点值表达式–>一般形式）

其实正变换的实现就是下列的矩阵相乘，反正我是没看出来
$$\left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& ...& ({\omega }_n^0{)}^{n-1}\\ & & & \\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& ...& ({\omega }_n^1{)}^{n-1}\\ & & & \\ .& .& ...& .\\ .& .& ...& .\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& ...& ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]\times \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ .\\ .\\ .\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}f(w_n^0)\\ f(w_n^1)\\ f(w_n^2)\\ .\\ .\\ .\\ f(w_n^{n-1})\end{array}\right]$$

---

• Add
  快速傅里叶逆变化$IFFT$，将点值表达式转换为系数表达式
  一般都是用的系数表达式

---

矩阵相乘的第i行第j列等于
求和第一个矩阵的第i行的每一个数和第二个矩阵的第j列的每一个数的乘积

我们记$V$就为（系数矩阵）上列的第一个矩阵，接下来再定义一个矩阵$D$，
$$D=\left[\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& ...& ({\omega }_n^{-0}{)}^{n-1}\\ & & & \\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& ...& ({\omega }_n^{-1}{)}^{n-1}\\ & & & \\ .& .& ...& .\\ .& .& ...& .\\ ({\omega }_n^{-n+1}{)}^0& ({\omega }_n^{-n+1}{)}^1& ...& ({\omega }_n^{-n+1}{)}^{n-1}\end{array}\right]$$

---

那么计算矩阵$D\ast V$的$i,j$项，分为两种情况
①：$i=j$，(任何数除零外的零次方都为1)
$$(D\ast V{)}_{i,j}=\sum _{k=0}^{n-1}{D}_{i,k}\ast V_{k,j}=\sum _{k=0}^{n-1}{\omega }_n^{-ik}\ast {\omega }_n^{jk}=\sum _{k=0}^{n-1}{\omega }_n^{(j-i)k}=\sum _{k=0}^{n-1}{\omega }_n^0=n$$
②：$i\ne j$
$$(D\ast V{)}_{i,j}=\sum _{k=0}^{n-1}{D}_{i,k}\ast V_{k,j}=\sum _{k=0}^{n-1}{\omega }_n^{-ik}\ast {\omega }_n^{jk}=\sum _{k=0}^{n-1}{\omega }_n^{(j-i)k}$$
$$=({\omega }_n^{j-i}{)}^0+({\omega }_n^{j-i}{)}^1+({\omega }_n^{j-i}{)}^2+...+({\omega }_n^{j-i}{)}^{n-1}$$
发现这个公式是一个以${\omega }_n^{j-1}$为公比的等比数列，


那么就可以转换为
$$\frac{({\omega }_n^{j-i}{)}^0-({\omega }_n^{j-i}{)}^1\ast ({\omega }_n^{j-i}{)}^{n-1}}{1-({\omega }_n^{j-i}{)}^1}=\frac{1-({\omega }_n^{j-i}{)}^n}{1-{\omega }_n^{j-i}}=\frac{0}{1-{\omega }_n^{j-i}}=0$$
单位复根的$n$次方$=0$，见上文单位复根定义
因为此公式的前提是$j\ne i$，所以分母一定不为$0$

故$j\ne i$时，$(D\ast V{)}_{i,j}=0$

---

用$D\ast V\ast f$去转换为点值表达式，去带最上面的这一板块的公式，你会惊讶地发现
$$\left[\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& ...& ({\omega }_n^{-0}{)}^{n-1}\\ & & & \\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& ...& ({\omega }_n^{-1}{)}^{n-1}\\ & & & \\ .& .& ...& .\\ .& .& ...& .\\ ({\omega }_n^{-n+1}{)}^0& ({\omega }_n^{-n+1}{)}^1& ...& ({\omega }_n^{-n+1}{)}^{n-1}\end{array}\right]\times \left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& ...& ({\omega }_n^0{)}^{n-1}\\ & & & \\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& ...& ({\omega }_n^1{)}^{n-1}\\ & & & \\ .& .& ...& .\\ .& .& ...& .\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& ...& ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]\times \left[\begin{array}{c}f(w_n^0)\\ f(w_n^1)\\ f(w_n^2)\\ .\\ .\\ .\\ f(w_n^{n-1})\end{array}\right]=$$
$$\left[\begin{array}{ccccc}n& 0& 0& ...& 0\\ 0& n& 0& ...& 0\\ 0& 0& n& ..& 0\\ ...& ...& ...& ...& 0\\ 0& 0& 0& ...& n\end{array}\right]\times \left[\begin{array}{c}f(w_n^0)\\ f(w_n^1)\\ f(w_n^2)\\ .\\ .\\ .\\ f(w_n^{n-1})\end{array}\right]=\left[\begin{array}{c}n\ast a_0\\ n\ast a_1\\ n\ast a_2\\ .\\ .\\ .\\ n\ast a_{n-1}\end{array}\right]$$
所以最后对答案全部$/n$就是点值表达式了，这也是为什么我们为这么定义$D,V$了

逆变换就相当于把正变换过程中的${\omega }_n^k$换成$w_n^{-k}$，之后结果除以n就可以了——摘自某dalao博客

离散傅里叶变换实现

理论

之前的思路全都是递归思想，实现出来后发现吓死个人，所以我们考虑转成迭代
以下的图摘自学长大佬：

学长让我们换成二进制看看：

可以发现终序列是原序列每个元素的翻转。
于是我们可以先把要变换的系数排在相邻位置，从下往上迭代。
在这里给出一个参考的方法：
我们对于每个 i，假设已知 i-1 的翻转为 j。考虑不进行翻转的二进制加法怎么进行：从最低位开始，找到第一个为 0 的二进制位，将它之前的 1 变为 0，将它自己变为 1。因此我们可以从 j 的最高位开始，倒过来进行这个过程
——摘自某dalao的博主

所以我们才会把这个$FFT$跟蝴蝶操作搞在一起，盗一波百度的图片

模板

void FFT ( complex *c, int f ) {
	for ( int i = 0;i < len;i ++ )
		if ( i < r[i] )
			swap ( c[i], c[r[i]] );
	for ( int i = 1;i < len;i <<= 1 ) {
		complex omega ( cos ( pi / i ), f * sin ( pi / i ) );
		for ( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w ( 1, 0 );
			for ( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = c[j + k], y = w * c[j + k + i];
				c[j + k] = x + y;
				c[i + j + k] = x - y;
			}
		}
	}
}

---

模板的板题运用

例题：洛谷P3803【模板】多项式乘法（FFT）

题目

递归版CODE

#include <cmath>
#include <cstdio>
using namespace std;
#define MAXN 3000005
struct complex {
	double real, i;
	complex () {}
	complex ( double xx, double yy ) {
		real = xx;
		i = yy;
	}
}a[MAXN], b[MAXN];

complex operator + ( complex s, complex t ) {
	return complex ( s.real + t.real, s.i + t.i );
}
complex operator - ( complex s, complex t ) {
	return complex ( s.real - t.real, s.i - t.i );
}
complex operator * ( complex s, complex t ) {
	return complex ( s.real * t.real - s.i * t.i, s.real * t.i + s.i * t.real );
}

const double pi = acos ( -1.0 );

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 w = complex ( cos ( 2 * pi / limit ), inv * sin ( 2 * pi / limit ) ), p = complex ( 1, 0 );
	for ( int i = 0;i < ( limit >> 1 );i ++, p = p * w ) {
		a[i] = a1[i] + p * a2[i];
		a[i + ( limit >> 1)] = a1[i] - p * a2[i];
	}
}

int main() {
	int n, m;
	scanf ( "%d %d", &n, &m );
	for ( int i = 0;i <= n;i ++ )
		scanf ( "%lf", &a[i].real );
	for ( int i = 0;i <= m;i ++ )
		scanf ( "%lf", &b[i].real );
	int limit = 1;
	while ( limit <= n + m )
		limit <<= 1; 
	FFT ( limit, a, 1 );
	FFT ( limit, b, 1 );
	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].real / limit + 0.5 ) );
	return 0;
}

迭代版CODE

推荐使用递推版，要比递归版快

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
#define maxn 3000005
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) {
		x = X, i = I;
	}
}A[maxn], B[maxn];
int len = 1;
int r[maxn];

double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

void FFT( complex *c, int f ) { //f=1系数转化为点值表达式w^1  f=-1点值转化为系数表达式w^(-1)
/*
蝴蝶发现：终序列是原序列每个元素二进制的翻转 
*/ 
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( c[i], c[r[i]] ); 
	for( int i = 1;i < len;i <<= 1 ) { //枚举迭代区间长度的一半 
		complex omega( cos( pi / i ), f * sin( pi / i ) );//区间长度本来是2i 就是要分成2i份 每一份是2pi/2i=pi/i 
		for( int j = 0;j < len;j += ( i << 1 ) ) {//枚举每一次迭代区间的开头
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
/*
只枚举迭代区间的左半部分
左半部分和右半部分进行计算
就可以算出上一层 直接覆盖即可
(w^k)^2=[w^(k+n/2)]^2 
左半部分是按照偶数分类
右半部分是按照奇数分类
f(x)=x*f1(x^2)+f2(x^2)
f1是奇数分类 f2是偶数分类 
*/ 
				complex x = c[j + k], y = w * c[j + k + i];
				c[j + k] = x + y;
				c[j + k + i] = x - y;
			}
		}
	}
}

int main() {
	int n, m;
	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 l = 0;
	while( len <= n + m ) {
		len <<= 1;
		l ++;
	}
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
/*
在原序列中i与i/2的关系是:i可以看做是i/2的二进制上的每一位左移一位得来
那么在反转后的数组中就需要右移一位
因为i直接左移一位
那么i二进制的右边第一位是没有考虑到的
那么如果那一位是1
反转后就应该是最高位为1 
*/
	FFT( A, 1 ); 
	FFT( B, 1 );
	for( int i = 0;i < len;i ++ )
		A[i] = A[i] * B[i];
	FFT( A, -1 );
	for( int i = 0;i <= n + m;i ++ )
		printf( "%d ", int( A[i].x / len + 0.5 ) );
	return 0;
}

给个版权吧：以上内容部分学习于

https://www.cnblogs.com/Tiw-Air-OAO/p/10162034.html
学校的lucky学长（没找到blog）
老师专讲
叉姐