fft详解

FFT详解


一些前置知识

任意角

详见高中必修4-1.1.1html

定义:按逆时针方向旋转造成的角叫作正角;按顺时针方向旋转造成的角叫作负角ios

弧度制

详见高中必修4-1.1.2数组

**定义:**把长度等于半径长的弧所对的圆心角叫1弧度的角,用符号 rad 表示(通常省略不写),读做弧度函数

常见弧度数:$360^\circ = 2\pi $ $180^\circ = \pi$优化

多项式

有$n$次多项式$F(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n$url

显然,它有$n+1$项spa

那么该多项式的系数表示就是$(a_0,a_1,a_2,\cdots,a_n)$.net

按照正常作法,在系数表示下,两个多项式相乘须要$O(n^2)$的时间复杂度code

用函数、方程的思想来看$F(x)$,至少肯定函数上$(n+1)$个点就能够推出惟一的$n$次函数$F(x)$orm

这$n+1$个点的坐标集合${(x_1,y_1),(x_2,y_2),\cdots,(x_{n+1},y_{n+1})}$称为该多项式的点值表示

在点值表示下,两个多项式相乘只需将对应y坐标相乘

又由于两个$n$项的多项式相乘所得的多项式有$2n+1$项

因此先将$A(x)$和$B(x)$的点值表达扩展到$2n+1$个点(上面说到的$n+1$个点是至少,实际上还能够再多)

再进行点值表达下的多项式相乘

以下:

$A(x)={(x_1,y_1),(x_2,y_2),\cdots,(x_{n+1},y_{n+1})}$

$B(x)={(x_1,y_1'),(x_2,y_2'),\cdots,(x_{n+1},y_{n+1}')}$

扩展后:

$A(x)={(x_1,y_1),(x_2,y_2),\cdots,(x_{2n+1},y_{2n+1})}$

$B(x)=(x_1,y'1),(x_2,y'2),\cdots(x{2n+1},y{2n+1})$

$C(x)=A(x)B(x)={(x_1,y_1y_1'),(x_2,y_2y_2'),\cdots,(x_{2n+1},y_{2n+1}y_{2n+1}')}$

可见,在点值表示下,两个多项式相乘只须要$O(n)$的时间复杂度

向量

详见高中必修4-2.1

**有向线段:**带有方向的线段。以$A$为起点,$B$为终点的有向线段记做$\vec{AB}$,起点写在终点前面;它的长度记做$|\vec{AB}|$

**三个要素:**起点、方向、长度。知道了有向线段的三个要素,它的终点就惟一肯定。

向量通常用有向线段表示

复数

定义:记$i=\sqrt{-1}$,把形如$a+bi$($a$、$b$均为实数)的数称为复数,其中$a$是实部,$b$是虚部

复平面:复数的平面中,$x$轴称为实轴,表示复数的实部;$y$轴称为虚轴,表示复数的虚部。为了方便理解,通常用从原点到点$(a,b)$的向量表示复数$a+bi$

模长:在复平面上,原点到点$(a,b)$的距离$\sqrt{a^2+b^2}$为复数$a+bi$的模长

辐角:在复平面上,实轴与表示复数$a+bi$的向量所造成的正角$\theta$称为复数$a+bi$的辐角

**运算:**复数运算与实数运算类似,以下: $$ (a+bi)\pm(c+di)=(a\pm c)+(b\pm d)i $$

$$ (a+bi)(c+di)=ac+adi+bci+bdi^2=ac+adi+bci+(-bd)=(ac-bd)+(ad+bc)i $$

$$ \frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(ac-bd)+(bc-ad)i}{c^2+d^2} $$

注意!!!(敲黑板*2)有个结论:复数乘法,辐角相加,模长相乘

欧拉公式

公式长这个样子:$e^{ix}=cos (x) + i sin(x)$。其中$x$是一个实数

这篇博客中有详细的证实

单位根

数学上,n次单位根是n次幂为1的复数。它们位于复平面的单位圆上,构成n边形的顶点,其中一个顶点是1。(From Baidu)

**单位圆:**在复平面上,以原点为圆心,单位长度为半径所做出的圆

**$n$次单位根:**在复平面上,以原点为起点,单位圆的$n$等分点为终点的向量所表示的$n$个复数中,辐角最小的向量所对应的复数$\omega_n$,称为$n$次单位根。将与$x$轴正半轴重合的那个向量记为$\omega_n^0$(那个上标其实是乘方),沿逆时针方向将剩余的向量顺序标记为$\omega_n^1,\omega_n^2,\cdots,\omega_{n-1}$。因此$\omega_n^n=\omega_n^0=1$。

将单位根和欧拉公式摆在一块儿,将上面提到的欧拉公式中的$x$换成$2\pi$(表示$n$次单位根$\omega_n^1$的向量与实轴所造成的夹角所对应的单位圆上的弧的长度),能够获得:$$e^{2\pi i}=cos(2\pi)+i sin(2\pi)=1=\omega_n^n$$

因此:$e^{\frac{2\pi i}{n}}=cos(\frac{2\pi}{n})+isin(\frac{2\pi}{n})=\omega_n$

咱们称此时的单位根$\omega_n$为主次单位根

那么其它的单位根就有:$\omega_n^k=cos(\frac{2\pi k}{n})+isin(\frac{2\pi k}{n})(0\leq k\leq n)$

消去引理:$\omega_{dn}^{dk}=\omega_n^k$ ($n$、$d$、$k$均为整数且$n\geq0,d>0,k\geq0$)。证实以下:$$\omega_{dn}^{dk}=(e^{\frac{2\pi i}{dn}})^{dk}=(e^{\frac{2\pi i}{n}})^k=\omega^k_n$$

**折半引理:**若是$n>0$且$n$为偶数,那么$n$个$n$次单位复数根的平方的集合就是$\frac{n}{2}$g个$\frac n 2$次单位复数根的集合。证实以下: $$ \begin{equation} \begin{aligned} & 根据消去引理,咱们有(\omega_n^k)^2=\omega_n^{2k}且(\omega_n^{n+\frac k 2})^2=\omega_n^{2k+n}=\omega_n^{2k}=(\omega_n^k)^2\ & 所以,\omega_n^k与\omega_n^{k+\frac n 2}平方相同\ & 从而能够推出\omega_n^k=-\omega_n^{k+\frac n 2}\ \end{aligned} \end{equation} $$ 求和引理:$\sum^{n-1}{j=0}(\omega^k_n)^j=0$,证实以下: $$ \sum^{n-1}{j=0}(\omega^k_n)^j=\frac{(\omega^k_n)^n-1}{\omega^k_n-1}=\frac{(\omega^n_n)^k-1}{\omega^k_n-1}=\frac{1^k-1}{\omega^k_n-1}=0 $$

矩阵

**定义:**矩阵是一个按照长方阵列排列的实数或负数集合,以下大小为$n\times m$的矩阵$A$,其中第$i$行第$j$列的数记做$a_{ij}$ $$ \begin{equation} \begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1m}\ a_{21} & a_{22} & \cdots & a_{2m}\ \vdots & \vdots & \ddots & \cdots\ a_{n1} & a_{n2} & \cdots & a_{nm}\ \end{array} \end{equation} $$

**矩阵乘法:**矩阵运算中比较经常使用的是矩阵乘法。矩阵乘法只有在第一个矩阵的行数和第二个矩阵的列数相同时才有意义,所以矩阵乘法不知足交换律。矩阵乘法法则为:有$p\times n$的矩阵$A$和$m\times p$的矩阵$B$,同时有矩阵$C=A\times B$,那么$C_{ij}=\sum_{k=1}^{p}a_{ki}\times b_{jk}$

**范德蒙德矩阵:**一个$n\times n$阶的范德蒙德矩阵$V$由包含$n$个数的数列$a_0,a_1,\cdots,a_{n-1}$决定,其中$v_{ij}=a_j^i(0\leq i,j<n)$,即: $$ \begin{equation} V(a_1,a_2,\cdots,a_n)= \left( \begin{array}{cccc} (a_0)^0 & (a_1)^0 & \cdots & (a_{n-1})^0\ (a_0)^1 & (a_1)^1 & \cdots & (a_{n-1})^1\ (a_0)^2 & (a_1)^2 & \cdots & (a_{n-1})^2\ \vdots & \vdots & \ddots & \vdots\ (a_0)^{n-1} & (a_1)^{n-1} & \cdots & (a_{n-1})^{n-1}\ \end{array} \right)

\left( \begin{array}{cccc} 1 & 1 & \cdots & 1\ a_0 & a_1 & \cdots & a_{n-1}\ a_0^2 & a_1^2 & \cdots & a_{n-1}^2\ \vdots & \vdots & \ddots & \vdots\ a_0^{n-1} & a_1^{n-1} & \cdots & a_{n-1}^{n-1} \end{array} \right) \end{equation} $$ 范德蒙德矩阵很特殊,它的逆矩阵是这样的: $$ \left( \begin{array}{cccc} 1 & 1 & \cdots & 1\ a_0 & a_1 & \cdots & a_{n-1}\ a_0^2 & a_1^2 & \cdots & a_{n-1}^2\ \vdots & \vdots & \ddots & \vdots\ a_0^{n-1} & a_1^{n-1} & \cdots & a_{n-1}^{n-1} \end{array} \right) ^{-1}=\frac{1}{n} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1\ a_0 & a_1 & \cdots & a_{n-1}\ a_0^2 & a_1^2 & \cdots & a_{n-1}^2\ \vdots & \vdots & \ddots & \vdots\ a_0^{n-1} & a_1^{n-1} & \cdots & a_{n-1}^{n-1} \end{array} \right) $$


如何加速多项式乘法

通常的多项式乘法题目都是给出两个多项式的系数表示,求这两个多项式相乘后所得的多项式的系数表示

对于这样的多项式乘法,最简单的方法就是两个for循环,也是系数表示下的多项式乘法,时间复杂度$O(n^2)$。代码以下,很是简洁直观:

for(int i=1;i<=n;i++){
	for(int j=1;j<=m;j++){
		ans[i+j-1]=a[i]*b[j];
	}
}

太慢了!!!

上面还提到,点值表示下的多项式乘法只须要$O(n)$

那若是咱们能够很快地将一个多项式的系数表示转为点值表示,而后在点值表示下进行多项式乘法,最后很快地将相乘获得的多项式点值表示转为系数表示,咱们就能够在多项式乘法中有一个优秀的时间复杂度了

这里介绍两种变换:

**DFT:**离散傅里叶变换,全称Discrete Fourier Transform,是指将多项式从系数表示转为点值表示的过程

**IDFT:**离散傅里叶逆变换,全称Inverse Discrete Fourier Transform,是指将多项式从点值表示转为系数表示的过程


DFT&FFT

DFT

将系数表示转为点值表示最简单的作法就是对多项式$A(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n$中的$x$取$n+1$个不一样的值,逐个代入多项式,求出对应的$y$

显然,将单个$x$代入求$y$须要$O(n)$的时间,那么将$n+1$个$x$代入求对应$n$个$y$就须要$O(n^2)$的时间复杂度

仍是太慢了。。。

这个时候就用到FFT加速DFT

FFT

**FFT:**快速离散傅里叶变换,全称Fast Fourier Transformation,就是加速后的DFT

在FFT中,采起分治的策略优化时间

为了方便分治,通常把多项式的项数补成2的正整数次幂,即将本来没有的更高次项的系数取0

下文默认$n$为2的正整数次幂

有$n$项多项式$A(x)=a_0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1}$

对这个多项式进行奇偶分治,获得如下两个多项式:

$F(x)=a_1x+a_3x^3+\cdots+a_{n-2}x^{n-2}$

$G(x)=a_0+a_2x^2+\cdots+a_{n-1}x^{n-1}$

每项的次数看上去挺不友善的。。。想办法化简一下,获得如下两个多项式:

$F'(x)=a_1+a_2x+a_5x^2+\cdots+a_{n-2}x^{\frac{n-3}2}$

$G'(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-1}x^{\frac{n-1}2}$

上面几个多项式的关系是这样的:

$F(x)=xF'(x^2)$

$G(x)=G'(X^2)$

$A(x)=F(x)+G(x)=xF'(x^2)+G'(x^2)$

而后为了方便FFT优化,选取的$n$个$x$为$n$次单位根$\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$

注意!!!是$n$个而不是$n+1$个,由于这里的多项式项数为$n$,即最高次项的次数只到$n-1$,因此只用$n$个点就能够肯定这个多项式

将$x=\omega_n^k$代入上式,获得:(消去引理) $$ \begin{equation} \begin{aligned} A(\omega_n^k)&=\omega_n^kF'(\omega_n^{2k})+G'(\omega_n^{2k})\ &=\omega_n^kF'(\omega_{n/2}^k)+G'(\omega_{n/2}^k) \end{aligned} \end{equation} $$ 由于在单位复数根$\omega_b^a$中,$a\leq b$

因此上式只能用于$k\leq n/2$的状况

设$k'\leq n/2$,那么当$k>n/2$时有$k=k'+n/2$,代入,获得: $$ \begin{equation} \begin{aligned} A(\omega_n^k)&=\omega_n^kF'(\omega_n^k)+G'(\omega_n^k)\ &=\omega_n^{k'+n/2}F'(\omega_n^{k'+n/2})+G'(\omega_n^{k'+n/2})\ &=\omega_n^{k'+n/2}F'(\omega_n^{2k'+n})+G'(\omega_n^{2k'+n})\ &=\omega_n^{k'+n/2}F'(\omega_n^{2k'}\times\omega_n^n)+G'(\omega_n^{2k'}\times\omega_n^n)\ &=\omega_n^{k'+n/2}F'(\omega_n^{2k'})+G'(\omega_n^{2k'})\ &=\omega_n^{k'+n/2}F'(\omega_{n/2}^{k'})+G'(\omega_{n/2}^{k'})\ &=-\omega_n^{k'}F'(\omega_{n/2}^{k'})+G'(\omega_{n/2}^{k'}) \end{aligned} \end{equation} $$ 由于$k$和$k'$都小于等于$n/2$,因此能够当作同一个东西

对比上面两个结果,就能够发现两个式子只有一个常数项不一样

因此每次只须要分治下去,求出$F'(\omega_{n/2}^k)$和$G'(\omega_{n/2}^k)$(注意这两个是不同的,它们每项的系数不同)

而后就能够合并为$A(\omega_n^k)$

即: $$ \left{ \begin{aligned} &a_{k} =\omega_n^kf'_k+g'k\ &a{k+n/2} =-\omega_n^kf'_k+g'_k \end{aligned} \right. $$ 这个分治的过程当中,每次合并须要$O(n)$,一共分治$logn$次,因此时间复杂度为$O(nlogn)$

更多的优化——迭代

上面的分治很容易想到用递归能够实现

可是递归很是耗内存,又跑得慢

因此尝试用迭代代替递归以让FFT更优秀

问题来了,怎样迭代呢?

找一下规律,就能够发现,奇偶分治到最底层($n=1$)后,结果序列的二进制形式为原序列的二进制翻转后结果(好像不是特别准确)

直观一点,举个例子:

当$n=8$时,进行奇偶分治:

(0 1 2 3 4 5 6 7)

(0 2 4 6)(1 3 5 7)

(0 4)(2 6)(1 5)(3 7)

(0)(4)(2)(6)(1)(5)(3)(7)

原序列,结果序列的二进制对比:

原序列 0 1 2 3 4 5 6 7
原序列二进制 000 001 010 011 100 101 110 111
结果序列 0 4 2 6 1 5 3 7
结果序列二进制 000 100 010 110 001 101 011 111

规律就很明显了,就是二进制反过来

具体代码实现中的二进制翻转推导能够看下这篇文,很好理解

因此就能够直接跳过从原序列分治到单项的过程,直接用分治完的结果序列合并上去

成功优化


IDFT&IFFT

若是把DFT看做矩阵乘法,能够获得: $$ \begin{equation} \left( \begin{array}{ccccc} (x_0)^0 & (x_0)^1 & (x_0)^2 & \cdots & (x_0)^{n-1}\ (x_1)^0 & (x_1)^1 & (x_1)^2 & \cdots & (x_1)^{n-1}\ (x_2)^0 & (x_2)^1 & (x_2)^2 & \cdots & (x_2)^{n-1}\ \vdots & \vdots & \vdots & \ddots & \vdots\ (x_{n-1})^0 & (x_{n-1})^1 & (x_{n-1})^2 & \cdots & (x_{n-1})^{n-1} \end{array} \right) \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \end{equation} $$ 将对应$x_0,x_1,x_2,\cdots,x_{n-1}$的$\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$带入,获得: $$ \begin{equation} \left( \begin{array}{ccccc} (\omega_n^0)^0 & (\omega_n^0)^1 & (\omega_n^0)^2 & \cdots & (\omega_n^0)^{n-1}\ (\omega_n^1)^0 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1)^{n-1}\ (\omega_n^2)^0 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2)^{n-1}\ \vdots & \vdots & \vdots & \ddots & \vdots\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1})^{n-1} \end{array} \right) \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \end{equation} $$ 化简,获得: $$ \begin{equation} \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1}\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{array} \right) \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \end{equation} $$

那么能够获得: $$ \begin{equation} \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right)

\left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{1} & \omega_n^{2} & \cdots & \omega_n^{(n-1)}\ \omega_n^0 & \omega_n^{2} & \omega_n^{4} & \cdots & \omega_n^{2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{(n-1)} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{array} \right) ^{-1} \end{equation} $$ 不难发现,须要求逆的矩阵是一个范德蒙德矩阵,那么它就能够这样转换一下: $$ \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{1} & \omega_n^{2} & \cdots & \omega_n^{(n-1)}\ \omega_n^0 & \omega_n^{2} & \omega_n^{4} & \cdots & \omega_n^{2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{(n-1)} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{array} \right) ^{-1} =\frac 1 n \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)}\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)} \end{array} \right) $$ 因此上面的式子就能够转化为: $$ \begin{equation} \left( \begin{array}{c} a_0\ a_1\ a_2\ \vdots\ a_{n-1} \end{array} \right) =\frac 1 n \left( \begin{array}{c} y_0\ y_1\ y_2\ \vdots\ y_{n-1} \end{array} \right) \left( \begin{array}{ccccc} \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n-1)}\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n-1)}\ \vdots & \vdots & \vdots & \ddots & \vdots\ \omega_n^0 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)(n-1)} \end{array} \right) \end{equation} $$ 因此只要对IDFT采起相似FFT的奇偶分治优化,就能将$O(n^2)$的IDFT加速为$O(nlogn)$的IFFT

具体为将FFT中的单位复数根$\omega_n^k$取为$\omega_n^{-k}$

因此能够将FFT和IFFT合到同一个函数中,区别为传入参数中的type为1仍是-1

完成(撒花


例题:luogu P3803 【模板】多项式乘法(FFT)

代码以下,详见注释:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

using namespace std;
const int N=3000000; 
const double pi=acos(-1.0);
int n,m,nn,bit,rev[N];

struct Complex{//定义复数类以及它的运算 
	double a,b;
	Complex(){}
	Complex(double aa,double bb){
		a=aa;b=bb;
	}
	Complex operator + (const Complex& x) const{
		return Complex(a+x.a,b+x.b);
	}
	Complex operator - (const Complex& x) const{
		return Complex(a-x.a,b-x.b);
	}
	Complex operator * (const Complex& x) const{
		return Complex(a*x.a-b*x.b,a*x.b+b*x.a);
	}
}a[N],b[N];

void init(){//初始化 
	for(int i=0;i<N;i++) a[i].a=a[i].b=b[i].a=b[i].b=0;
	nn=1;
	bit=0;
}

void before_fft(){
	while(nn<=n+m+1){ //将多项式项数扩展到2的正整数次幂 
		nn<<=1;
		bit++;
	}
	for(int i=0;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//求对应二进制翻转 
}

void fft(Complex* a,int typ){
	for(int i=0;i<nn;i++){//将数组变到奇偶分治到底层时的顺序 
		if(i<rev[i]) swap(a[i],a[rev[i]]); //if不能少,不然会有数交换两次后位置不变 
	}
	for(int i=2;i<=nn;i<<=1){//每次要合并的两个多项式的项数和 
		Complex wn=Complex(cos(2*pi/i),typ*sin(2*pi/i));//经过欧拉公式求主次单位根 
		for(int st=0;st<nn;st+=i){//枚举归并位置 
			Complex wnk=Complex(1,0); 
			for(int k=0;k<i/2;k++){//枚举w_n的幂数k 
				Complex x=a[st+k];//对应fft推导过程当中的g'k 
				Complex y=wnk*a[st+k+i/2];//对应fft推导过程当中的wnk f'k
				a[st+k]=x+y;//对应a_k 
				a[st+k+i/2]=x-y;//对应a_k+n/2 
				wnk=wnk*wn;
			}
		}
	} 
	if(typ==-1){
		for(int i=0;i<=nn;i++) a[i].a/=nn;//根据ifft的推导,这里要记得除以nn 
	}
}

int main(){
	init();
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&a[i].a);
	for(int i=0;i<=m;i++) scanf("%lf",&b[i].a);
	before_fft();
	fft(a,1);
	fft(b,1);
	//fft后的a,b分别为这两个多项式在点值表达下当x取wn1,wn2,...,wn(n-1)时的y 
	for(int i=0;i<=nn;i++) a[i]=a[i]*b[i];//点值表达下的多项式乘法,只将对应y相乘 
	fft(a,-1);//ifft 
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].a+0.5));//这里要四舍五入 
	printf("\n");
	return 0;
}
相关文章
相关标签/搜索