FFT学习笔记

FFT(Fast Fourier Transform)

快速傅里叶变换

---

引入

• 百度一下

这里的很多的工程上的文章和学术性的文章都是以音频和图像的处理和生活中的实际用途还有纯数学角度讲的，但是作为一个OIer，我们一般是不会用到这些的，所以这里就不讲解什么时域频域转换和波的分析之类的东西，而是将如何在OI中应用。

概念与前置

人物了解

• 傅里叶

其实这里要讲的傅里叶变换和完整的傅里叶变换有一定的区别，我们下面来看看傅里叶变换的两个公式:

• 傅里叶变换
  $$F(\omega )=\mathcal{F}[f(t)]=\underset{-\infty }{\overset{\infty }{\int }}f(t)e^{-i\omega t}dt$$

• 傅里叶逆变换
  $$f(t)={\mathcal{F}}^{-1}[F(\omega )]=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}F(\omega )e^{i\omega t}d\omega$$

其中的$\mathcal{F}$代表傅里叶变换，${\mathcal{F}}^{-1}$表示傅里叶逆变换。

下面扯一些和算法并无多大联系的知识。

如果不想看直接可以调到下一个标题

这两个变换公式中，$f(t)$为以$t$为周期的周期函数，且$t$还满足狄利赫里条件。

狄利赫里条件
内容概括:也就是在给定的周期函数在有限的区间内,只有有限个第一类间断点和有限个极大值和极小值。满足该条件就符合傅里叶级数收敛。

感觉每次学啥都是递归似的学习啊QWQ

傅里叶级数
内容概括:任何周期函数都可以用正弦函数和余弦函数构成的无穷级数来表示（选择正弦函数与余弦函数作为基函数是因为它们是正交的），后世称傅里叶级数为一种特殊的三角级数，根据欧拉公式，三角函数又能化成指数形式，也称傅立叶级数为一种指数级数。

• 傅里叶级数的收敛

设$f$是周期为$2\pi$且在$[-\pi ,\pi ]$上的可积或绝对可积的函数，如果$f$在$x_0$处存在倒数$f^{{}^{\prime }}(x_0)$，或是有两个有限的单侧倒数:
$$f_{+}^{{}^{\prime }}(x_0)=\underset{t\to 0^{+}}{\lim }\frac{f(x_0+t)-f(x_0)}{t}$$
$$f_{-}^{{}^{\prime }}(x_0)=\underset{t\to 0^{+}}{\lim }\frac{f(x_0-t)-f(x_0)}{-t}$$
那么$f$的$\mathrm{F}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{r}$级数在$x_0$处收敛于$f(x_0)$，如果$f$在$x_0$处仅有两个有限的广义单侧导数:
$$\underset{t\to 0^{+}}{\lim }\frac{f(x_o+t)-f(x_0+0)}{t}$$
$$\underset{t\to 0^{+}}{\lim }\frac{f(x_0-t)-f(x_0-0)}{-t}$$
那么$f$的$\mathrm{F}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{r}$级数在$x_0$处收敛于$\frac{1}{2}(f(x_0+0)+f(x_0-0))$.

参考文章

傅里叶级数的三角式:
$$f(t)=a_0+\sum _{n=1}^{\infty }(a_n\cos (n\omega t)+b_n\sin (n\omega t))$$

其中$\omega =\frac{2\pi }{T}$，$T$为函数周期。$a_n$和$b_n$分别控制了正弦波的振幅与频率。说好的不讲这些的呢？

它还可以用复指数形式和积分形式来表示:
$$f(t)=\sum _{n=-\infty }^{\infty }{F}_n{e}^{in\omega t}$$
$$F_n=\frac{1}{T}{\int }_0^{T}f(t)e^{-in\omega t}dt$$
其中的$F$就是周期函数$f$的傅里叶级数(Fourier Series, FS)

---

再回到傅里叶变换的两个式子，里面的$F(\omega )$称为$f(t)$的象函数，而$f(t)$被称为$F(\omega )$的原象函数。

象函数与原象函数
拉普拉斯变换
内容概述:这是一种线性的积分变换，可以将一个有参数实数$t(t\ge 0)$的函数变换成一个参数为复数的函数$S$。
在拉普拉斯变换中对于$f(t)$的双边拉普拉斯变换$F(s)$,我们称其中的$F(\omega )$为$f(t)$的象函数，而$f(t)$则为$F(\omega )$的原象函数。

---

跳出递归

下面开始由最简单的入手一步步讲解$\mathrm{F}\mathrm{F}\mathrm{T}$

多项式

• 定义

这个应该来看的人都知道吧

我们把形如$a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$的式子称作多项式(其中$a_0$可以形式化的写作$a_0{x}^0$)。
其中的$a_0$~$a_1$叫做多项式的系数，我们将一个关于$x$的多项式可以简记为$A(x)$。

多项式的通式表达:
$$A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$

一些规定:

• 我们令$n$为一个多项式的次数界，即次数范围，那么这个多项式中的任意一项的次数都要小于$n$。

• 我们定义两个多项式的加减法如下(两个多项式加起来的次数界去较大的一个):
  $$A(x)\pm B(x)=\sum _{i=0}^{n-1}(a_i\pm b_i)x^i$$

这个可以$O(n)$的简单计算啦。

• 我们接下来定义两个多项式的卷积(好像也就是乘法吧):

$$A(x)⨂B(x)=\sum _{i=0}^{A_n-1}\sum _{j=0}^{B_n-1}{a}_i{b}_i\cdot x^{i\times j}$$

两个多项式的卷积的次数界为$(A_n+B_n)$。

卷积式子也可以写成以下形式:
$$A(x)⨂B(x)=\sum _{i=0}^{n-1}\left(\sum _{j=0}^i{a}_j{b}_{i-j}\right)x^i$$

那么从定义式的角度来看，朴素的算法只能$O(n^2)$的计算，那么有没有更快一点的计算呢?我们接着看。

---

多项式的表达方法

• 第一种:系数表达法

也就是我们上面的定义中的方法，用一组系数$a_0$~$a_1$来表示一个多项式，我们可以将这个系数看作一个$n-1$维的向量，记作$a^{\to }=(a_0,a_1,\cdots ,a_n-1)$ ,然后向量一般可以与矩阵联系起来(后面有用)。

• 第二种:点值表达法

因为一个多项式$A(x)$的实质是关于$x$的一个函数，而一个$n-1$次的函数可以由$n$个点对$(x,y)$确定其形式形态，所以我们可以运用这个性质来表示一个多项式。对于一个次数界为$n-1$的多项式，我们可以选取$n$个$x$，记作$(x_0,x_1,\cdots ,x_{n-1})$，然后带入这个多项式，可以对应求出$n$个$y$，我们把它记作$(y_0,y_1,\cdots ,y_{n-1})$，那么现在我们就得到了这个函数(多项式)上的$n$个点对，我们就可以用这$n$个点对$(x_i,y_i)$来表示这个多项式。其实这也是离散化的表示一个函数(多项式)的方法。

然后我们将其形式化:

令 α = { ( x i , y i ) [ y i = A ( x i ) ] ∣ i ∈ [ 0 , n − 1 ] , i ∈ Z } \alpha=\{(x_i,y_i)[y_i=A(x_i)]|i\in [0,n-1],i\in Z\} α={(xi​,yi​)[yi​=A(xi​)]∣i∈[0,n−1],i∈Z}
这个$\alpha$即为$A(x)$的点值表达式。

我们用这个点值表达式有什么用处呢?
点值表达式的方便之处在于我们将两个多项式相乘如果用系数表达式则是$O(n^2)$的，而使用点值表达式，则只用把两个多项式对应$x_i$的点相乘即可，即变为$O(n)$的了。(这里的点相乘会用到复数的知识)。

复数

复数其实就是一种数域的扩展，为了解决一些实数解决不了的和平面几何等的问题。

• 百度一下

• 有关复数的运算

于是上面的点便可以定义为复数域$C$上的复数。
那么对于一个点$p$，我们可以将其表示为$p=x+yi$，其中$i$就是复数单位，其中$i$满足$i^2=-1$。

那么这个点的四则运算就可以定义如下(结合复数的运算):

对于两个点$p,q$，其中$p=x_p+y_{p}i,q=x_q+y_{q}i$，我们可以用$(x_p,y_p)$来表示点$p$。

• 加法$+$
  $$p+q=(x_p+x_q,y_p+y_q)$$

• 减法$-$
  $$p-q=(x_p-x_q,y_p-y_q)$$

• 乘法$\times$
  $$p\times q=(x_p\times x_q-y_p\times y_q,x_p\times y_q+x_q\times y_p)$$

• 除法$÷$
  $$\frac{p}{q}=\left(\frac{x_p\times x_q+y_p\times y_q}{x_q^2+y_q^2},\frac{y_p\times x_q-y_q\times x_p}{x_q^2+y_q^2}\right)$$

• 共轭复数$conj$【百度百科】
  $conj(p)=(x_p,-y_p)$

复数就介绍到这，我们的点值表达式中的点的运算就是如此（把点用复数表示了（复平面））。

---

多项式求值

那么我们对于一个$x_i$如何快速求出$A(x)$呢?
会的大佬可以跳过 Orz

先观察定义式（次数界为$n$）:
$$A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$

最笨的方法$O(n^2)$，每次暴力计算$x^i$。
优化一点就是$O(nlo{g}_{2}n)$,利用快速幂。
再好一点的就是$O(n)$预处理$x^i$,然后$O(n)$计算。

那么有没有极致一点的直接$O(n)$的方法呢？

我们就是用秦九韶算法(Horner法则)。
其实他的本质类似于分治递归处理。

我们将这个多项式展开:
$$A(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1}$$
每次提取一个$x$，然后一直这么做下去，就得到如下变换。
$$a_0+x\cdot (a_1+a_{2}x+a_3{x}^2+\cdots +a_{n-1}{x}^{n-2})$$
$$a_0+x\cdot (a_1+x\cdot (a_2+a_{3}x+\cdots +a_{n-1}{x}^{n-3}))$$
$$\cdots$$
的到最里面一层就是
$$\cdots (a_{n-2}+a_{n-1}x)$$
我们发现变成了一个一次的式子，那么将$x_i$带入，求出值为$v_0$，继续得到如下式子:
$$\cdots (a_{n-3}+v_{0}x)$$
然后同样计算，直到最后一个，那么每次都是一个一次的式子，$O(1)$的计算，然后总的复杂度就是$O(n)$啦。

代码实现就很简单啦

int v=a[n-1];
for(int i=n-2;i>=0;i--){
	v=v*x+a[i];
}

$smalltrick$

---

接下来我们就继续看多项式卷积，我们有如下卷积式子:
$$C(x)=A(x)⨂B(x)=\sum _{i=0}^{A_n-1}\sum _{j=0}^{B_n-1}{a}_i{b}_i\cdot x^{i\times j}$$

表示为点值表达式就有 C ( x ) = { x i , A ( x i ) ⋅ B ( x i ) ∣ i ∈ [ 0 , n − 1 ] , i ∈ Z } C(x)=\{x_i,A(x_i)\cdot B(x_i)|i\in [0,n-1],i\in Z\} C(x)={xi​,A(xi​)⋅B(xi​)∣i∈[0,n−1],i∈Z}

我们发现这样可以$O(n)$的计算(前面的点值表达方式说过)。

但虽然计算过程中复杂度十分优秀，但是预处理呢?

我们一般是得到一个系数表达式，然后转换为点值表达式的朴素复杂度为$O(n^2)$，先枚举$n$个$x_i$，然后每次$O(n)$的计算多项式的值。然后我们$O(n)$的运算，但大多数时候最终的答案要为系数表达式，所以最朴素的转变回来的方式就是$O(n^3)$高斯消元(也就是解一个知道$n$个点的$n$次方程)。

这个，，，，还不如暴力呢，emmmm。

那么有没有一种快速的变换来解决这种预处理，当然有啦，那就是-------------------------------> FFT

---

继续前置$\cdots \cdots$

插值-点值转换

我们继续定义一种对于多项式的运算，名叫$DFT(A(x))$（其实有的读者就应该知道这就是 离散傅里叶变换$discrete\ Fourier\ transform $）

它的作用就是将系数表达式的$A(x)$转换为点值表达式。

而我们同样的还要定义它的逆运算，为$IDFT(A(x))$，作用为将一个点值表达式转换为系数表达式的转换。

那么多项式的卷积就可以如下计算：

• $DFT(A(x))$ $DFT(B(x))$

• $C(x)={\sum }_{i=0}^{n-1}A(x_i)\cdot B(x_i)$

• $IDFT(C(x))$

最后的答案就为$C(x)$了。

但是如何快速进行$DFT$与$IDFT$呢？肯定不能用朴素的方法。

所以这里就要运用复数的神奇性质啦!

---

单位根

一种十分神奇且超级特殊的一类方程的根。

定义: 即对于满足方程$x^n=1$的$x$，被称为$n$次单位根。

但从表面来看，$n$为奇数时$x$只有1，偶数时有$-1,1$，这样的根太少，且用处不大，那么我们就将目光投向复数。

重要的性质 : $i^2=-1$,$i^3=-i$,$i^4=1,\cdots$

那么在复数上，$x^n=1$就几乎有无数多组解，那么这些解就被称为$n$次单位复数根。

因为有着$x^n=1$的性质，那么将其带入多项式，所有的$x^i$都可以变为1，那不就好解多啦。

如何求得复数意义下$x^n=1$的解呢？

• 欧拉公式

$e^{xi}=\cos x+i\sin x $

其中一个特殊的性质：$(\cos x+i\sin x{)}^k=\cos kx+i\sin kx$
证明：用三角恒等式，变形，然后推出$k$在奇数和偶数时的通式证明即可。

那么欧拉公式如何得到的呢？

证明：
通过$Taylor$展开得到（有高数基础的应该知道，不知道的也可以百度）：

$$e^x=1+x+\frac{x^2}{2!}+\cdots +\frac{x^n}{n!}+\cdots$$
$$\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots$$
$$\sin x=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots$$

然后我们将其扩展到复数域上，由于$i=\sqrt{-1}$

所以可以得到如下式子：
$$e^{xi}=1+xi-\frac{x^2}{2!}-\frac{x^3}{3!}i+\frac{x^4}{4!}+\frac{x^5}{5!}i-\cdots$$
变形得到：
$$=\left(1-\frac{x^2}{2!}+\frac{x^4}{4!}-\cdots \right)+i\cdot \left(x-\frac{x^3}{3!}+\frac{x^5}{5!}-\cdots \right)$$
$$=\cos x+i\sin x$$
于是得到证明。

如果我们将$x$用$\pi$来代替，就可以得到
$$e^{i\pi }=\cos \pi +i\sin \pi =-1+i\times 0=-1$$
即一个神奇的公式就是$e^{i\pi }+1=0$

然后我们接着将$x=2\cdot k\pi (k\in Z)$
那么又能够得到$e^{2\cdot k\pi i}=1$，那么这里的$e^{2\cdot k\pi i}$不就正好符合我们的方程$x^n=1$吗？
然后我们令$x=e^{\frac{2k\pi i}{n}}(k\in Z)$，然后我们记${\omega }_n^k=x$,那么这个就是一个$n$次单位复数根。
而且可以证明这样的$n$次单位复数根一定有$n$个：
证明：
令$k=nr+p,p\in [0,n-1]$，则${\omega }_n^k=e^{\frac{2k\pi i}{n}}=e^{\frac{2(nr+p)\pi i}{n}}=e^{\frac{2nr\pi i}{n}+\frac{2p\pi i}{n}}=e^{\frac{2nr\pi i}{n}}\times e^{\frac{2p\pi i}{n}}=1^r\times e^{\frac{2p\pi i}{n}}=e^{\frac{2p\pi i}{n}}$
有因为$p$有$n$个，那么这个根也就有$n$个。

记这些根为${\omega }_n^0$~${\omega }_n^{n-1}$,且它们在乘法意义下构成一个群，即满足${\omega }_n^i{\omega }_n^j={\omega }_n^{i+j}={\omega }_n^{(i+j)\mathrm{m}\mathrm{o}\mathrm{d}n}$。

还有一个特殊的性质是${\omega }_n^{-1}={\omega }_n^{n-1}$
证明:
$${\omega }_n^{-1}=1\times {\omega }_n^{-1}=e^{2\pi i}\times e^{\frac{-2\pi i}{n}}=e^{2\pi i-\frac{2\pi i}{n}}=e^{\frac{2(n-1)\pi i}{n}}={\omega }_n^{n-1}$$

---

三个引理

• 消去引理
  当$d̸=0$时，有${\omega }_{dn}^{dk}={\omega }_n^k$
  证明:${\omega }_{dn}^{dk}=e^{\frac{2dk\pi i}{dn}}=e^{\frac{2k\pi i}{n}}={\omega }_n^k$

• 折半引理
  当$n>0$且$n$为偶数，$({\omega }_n^k{)}^2={\omega }_{\frac{n}{2}}^k$,用文字描述就是$n$个$n$次单位根的平方的集合就等于$\frac{n}{2}$个$\frac{n}{2}$次单位根的集合。
  证明:$({\omega }_n^k{)}^2=(e^{\frac{2k\pi i}{n}}{)}^2=e^{\frac{2k\pi i}{n}\times 2}=e^{\frac{2k\pi i}{\frac{n}{2}}}={\omega }_{\frac{n}{2}}^k$

• 求和引理
  对于任意的$k\in [0,n]$，都有${\sum }_{j=0}^{n-1}({\omega }_n^k{)}^j=0$
  证明：用等比数列的和的公式，带入推导化简一下就出来了(公式太难打了就读者自己推一推吧QAQ)。

特殊的：${\omega }_n^n=1$

---

那么运用复数和单位根，可以将$DFT$与$IDFT$降到$O(nlo{g}_{2}n)$的复杂度，下面会继续详解方法。

---

马上接近正题了，加油！

蝴蝶变换

其实也叫$Cooley-Tukey$算法，它是将$radix-r$算法选取$2^k$为基，(其实基可以在数据的大小内任意选，但是一般选$2$来变换)，然后来做变换，也就是按照下标的$2^k$进制来进行分治，它一般具有$lo{g}_{2^k}n$的复杂度，但是只能处理大小为$2^k$的整倍数的数据。

例如我们以$2$为基，来将一个多项式$A(x)$的系数分成$a,b$两个集合。
令集合大小为$n,S=\frac{n}{2}$,那么就有： a = { 2 S a 1 + 2 S − 2 a 2 + ⋯ + 2 a S − 1 + a S } , b = { b 1 + 2 b 2 + ⋯ + 2 S b S } a=\{2^{S}a_1+2^{S-2}a_2+\cdots +2a_{S-1}+a_S\},b=\{b_1+2b_2+\cdots +2^{S}b_S\} a={2Sa1​+2S−2a2​+⋯+2aS−1​+aS​},b={b1​+2b2​+⋯+2SbS​}（其中的系数只是表示下标被分治的标志）。

那么就可以使用蝴蝶变换来神奇的实现该操作，具体形式如下图：

其实我这个图画的不好，（mspaint好难用QAQ）。

可以看这，详细的手动模拟

详细的学术论文

我在这就不多说了，蝴蝶逆变换也就是将${\omega }_n^k$改为$\frac{1}{2}{\omega }_n^{-k}$就是逆变换了。

逆变换的证明可以看这里超级好的大佬文章，用矩阵的方法证明。

那么代码上也就简单啦！

---

离散傅里叶变换

其实还有连续的傅里叶变换

利用前面的$cooley-Tukey$算法，我们一般将基设置为$2$，那么就是按照多项式的下标的奇偶性分类（前面提出过），所以我们记原多项式为$A(x)$，它被分成两个$\left(A_0(x),A_1(x)\right)$（偶，奇），我们可以将其写成如下形式:

$$A_0(x)=a_0+a_{2}x+a_4{x}^2+\cdots =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2\cdot i}(x^i)$$
$$A_1(x)=a_1+a_{3}x+a_5{x}^2+\cdots =\sum _{i=0}^{\frac{n}{2}-1}{a}_{2\cdot i+1}(x^i)$$
然后根据折半引理
$$A(x)=A_0(x^2)+x\cdot A_1(x^2)$$
我们将单位根带入，可知所有的$x^i$，那么直接用系数与单位根就可以计算了，并且根据折半引理，可以每次按奇偶折半分治下去，复杂度为$O(logn)$，每次变换$O(n)$，所以总的时间复杂度为$O(nlogn)$。

所以对于$k<\frac{n}{2}$，我们有：

$$A({\omega }_n^k)=A_0(({\omega }_n^k{)}^2)+{\omega }_n^k{A}_1(({\omega }_n^k{)}^2)$$
$$=A_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_1({\omega }_{\frac{n}{2}}^k)$$
$$=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}({\omega }_{\frac{n}{2}}^{ki})+{\omega }_n^k\sum _{i=0}^{\frac{n}{2}-1}{a}_{[2i+1]}({\omega }_{\frac{n}{2}}^{ki})$$

还有一个性质就是
$$A({\omega }_n^{\frac{n}{2}+k})=A_0({\omega }_{\frac{n}{2}}^k)+{\omega }_n^{\frac{n}{2}+k}{A}_1({\omega }_{\frac{n}{2}}^k)$$
又因为${\omega }_n^{\frac{n}{2}}=-1$
所以原式就为
$$A({\omega }_n^{\frac{n}{2}+k})=A_0({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_1({\omega }_{\frac{n}{2}}^k)$$
所以另一半也就可以求出了。

那么这个问题就较好的解决了，在$O(nlo{g}_{2}n)$的时间内。

---

离散傅里叶逆变换

$IDFT(InverseDiscreteFourierTransform)$

就是$FFT$的最后一步了(上一步转换为点值表达式后还要O(n)的对于项乘起来，然后在对于乘起来后的点值表达式再进行IDFT)。

IDFT它的实质就是解一个线性方程组，但不会像高斯消元一样，具体的说明可以看前面推荐过的这篇文章，因为数学公式太难打了

---

具体实现

根据蝴蝶变换和折半引理，分治递归处理是最直接的，但是有时效率较差，那么我们深入探究蝴蝶变换的过程发现，可以用迭代的方法，而基为$2$的变换，实质就是基于二进制位，但是数组顺序会在变换后打乱，所以我们要提前翻转一下，下面看代码，代码中有些解释。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
using namespace std;
const int M=2e6+1e5+10;
const db pi=acos(-1);
struct complex{
	db x,y;
	complex(db a=0,db b=0):x(a),y(b){}
	complex operator +(complex a){return complex(x+a.x,y+a.y);}
	complex operator -(complex a){return complex(x-a.x,y-a.y);}
	complex operator *(complex a){return complex(x*a.x-y*a.y,x*a.y+y*a.x);}
	//定义复数类，也可使用自带的complex，但是速度会比手写的慢。
}a[M],b[M],w[M];
int r[M],len,lg;

void DFT(complex *a,int n,int f){
//	for(int i=1,j=(n>>1),k;i<=n-2;i++){
//		if(i<j)swap(a[i],a[j]);
//		for(k=(n>>1);j>=k;j-=k,k>>=1);
//		if(j<k)j+=k;
//	}每次重新变换
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);//预处理变换后直接交换
	for(int i=2;i<=n;i<<=1){
		int now=i>>1;
		complex wn=w[i];wn.y*=f;
		//预处理了单位根，如果为逆变换，f=-1，然后将sin项乘以-1
		//不预处理的话精度误差会大一点，每次算的话，wn=complex(cos(2*pi/i),f*sin(2*pi/i))
		for(int j=0;j<n;j+=i){
			complex w=complex(1,0),x,y;
			for(int k=j;k<j+now;k++,w=w*wn){
				x=a[k],y=w*a[k+now];//枚举奇数项和偶数项，然后奇数项乘以单位根的次方
				a[k]=x+y;a[k+now]=x-y;//根据前面推导的式子进行加减
			}
		}
	}
	if(f==-1)for(int i=0;i<n;i++)a[i].x/=n;//其实逆变换每次除以2，可以简单写为最后除以n=2^s
}

void FFT(complex *a,complex *b,int n,int m){
	DFT(a,len,1);DFT(b,len,1);//先将两个多项式转换为点值表达式
	for(int i=0;i<=len;i++) a[i]=a[i]*b[i];//对应项相乘
	DFT(a,len,-1);//逆变换回来
	for(int i=0;i<=n+m;i++) printf("%d ",int(a[i].x+0.5));
	//输出系数，要注意精度误差
}
int x,n,m;
int main(){
//输入两个0~n和0~m的多项式系数，然后求其卷积后输出
//洛谷FFT模板
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%d",&x),a[i].x=x;
	for(int i=0;i<=m;i++) scanf("%d",&x),b[i].x=x;
	for(len=1;len<=n+m;len<<=1,++lg);
	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(lg-1));//预处理变换
	for(int i=2;i<=len;i<<=1)w[i]=complex(cos(2*pi/i),sin(2*pi/i));//预处理单位根
	FFT(a,b,n,m);
	return 0;
}

---

结束与例题

类似FFT的还有：

• NTT（快速数论变换）：可以解决FFT的精度问题，用于在取模意义下的卷积。

• FWT（快速沃尔什变换）：解决集合卷积的快速算法，【我的讲解】

• FMT（快速莫比乌斯变换）：和莫比乌斯反演有点关系，我也不清楚，好像也可以解决一些奇怪的集合卷积问题。

---

一些题目：

• 模板题

• 高精度乘法

• [ZJOI2014]力

• 残缺的字符串

• 万径人踪灭

• 九连环(高精快速幂)

• HDU4609 3-idiots

• 我自己的题抛色子（完全背包计数FFT优化）

---

如果有错的地方请指出，Thanks♪(･ω･)ﾉ。