FFT快速傅里叶变换详解

介绍

简而言之，这个东西用来做多项式乘法。

比如说，有两个多项式：
3x2+2x+1,2x2+x+43x^2+2x+1~,~2x^2+x+4 3x2+2x+1 , 2x2+x+4

那么他们乘起来就是
6x4+7x3+16x2+9x+46x^4+7x^3+16x^2+9x+4 6x4+7x3+16x2+9x+4

FFTFFTFFT 也就是做了这样一件事。

但是如果暴力取做这个乘法的话，时间复杂度是 O(n2)O(n^2)O(n2) 的，然而优秀的 FFTFFTFFT 只需要 O(nlogn)O(nlogn)O(nlogn)。

为了方便，下面设相乘的多项式为 A,BA,BA,B，乘起来得到的多项式为 CCC，也就是 A×B=CA \times B = CA×B=C，以及设 AAA 是一个 nnn 次函数，BBB 是一个 mmm 次函数，CCC 是一个 kkk 次函数。

（有什么不严谨的地方还请见谅qwq）

前置芝士

三角函数总得要会呀
得知道复数是个什么东西
平面直角坐标系得会用吧

平面直角坐标系和三角函数就不讲了。要是这个东西都不知道学什么 FFTFFTFFT……

复数

前置芝士——虚数： 定义虚数单位 i=−1i=\sqrt{-1}i=−1，虚数就是形如 a+bi(b≠0)a+bi~~(b\neq 0)a+bi (b=0) 的数，其中 a,ba,ba,b 都是实数。

复数： 形如 a+bia+bia+bi 的数，其中 a,ba,ba,b 都是实数。说白了，复数就是实数与虚数的并集。

实部： 也就是 aaa。

虚部： 也就是 bibibi。

复数加法： 实部与实部相加，虚部与虚部相加。例： (2+3i)+(−9+4i)=−7+7i(2+3i)+(-9+4i)=-7+7i(2+3i)+(−9+4i)=−7+7i。

复数减法： 与加法类似，实部与实部相减，虚部与虚部相减。

复数乘法：
根据多项式乘法的规则去乘即可。例： (2−2i)×(3+1i)=2×3+2i−6i−2i2=6−4i+2=8−4i(2-2i)\times(3+1i)=2\times3+2i-6i-2i^2=6-4i+2=8-4i(2−2i)×(3+1i)=2×3+2i−6i−2i2=6−4i+2=8−4i。别忘了 i=−1i=\sqrt{-1}i=−1，即 i2=−1i^2=-1i2=−1。
由此可以得出，复数乘法的规则是：实部与实部相乘 +++ 虚部与虚部相乘 === 积的实部；实部与虚部相乘 +++ 虚部与实部相乘 === 积的虚部。

根据这种计算规则，我们可以写出这样的复数结构体：

struct comp{double x,y;comp(double xx=0,double yy=0):x(xx),y(yy){}comp operator +(const comp b){return (comp){x+b.x,y+b.y};}comp operator -(const comp b){return (comp){x-b.x,y-b.y};}comp operator *(const comp b){return (comp){x*b.x-y*b.y,x*b.y+y*b.x};}
};

正题

我们发现，将上面两个多项式乘起来之后，得到了另一个多项式。假如我们将这三个多项式看成三个函数，那么会怎么样呢？

初中时就学过：n+1n+1n+1 个点可以确定一个 nnn 次函数。那我们只需要找出 k+1k+1k+1 个在函数 CCC 上面的点，那么就可以求出 CCC 了。

但是，直接找这些点的过程就是 O(n2)O(n^2)O(n2) 的了……

于是，FFTFFTFFT 对这个方面进行了优化！他不找实数点，而是去找复数点。

对于一个复数，我们可以将其表示为 a+bia+bia+bi，那么可以在平面直角坐标系中将其表示成一个坐标为 (a,b)(a,b)(a,b) 的点。

比如说，对于一个复数 3+4i3+4i3+4i，在平面直角坐标系中就是：

这个点有另外一个表示方法：我们设它与原点之间的线段的长度为 zzz，与 xxx 轴之间的夹角为 θ\thetaθ，也就是这样子：

因为 x=cos⁡(θ)×z,y=sin⁡(θ)×zx=\cos(\theta)\times z,y=\sin(\theta) \times zx=cos(θ)×z,y=sin(θ)×z，所以这个点的坐标又可以表示为：(cos⁡(θ)×z,sin⁡(θ)×z)(\cos(\theta) \times z,\sin(\theta) \times z)(cos(θ)×z,sin(θ)×z)。

单位圆

单位圆： 以平面直角坐标系的原点为圆心做一个半径为 111 的圆，这个圆我们称它为单位圆。

就是这样的一个东西：

单位根

定义： 我们称满足 xn=1x^n=1xn=1 的复数 xxx 为 nnn 次单位根。

当然，对于一个 nnn，是可以有多个单位根的，其中一个就是 111。

怎么找呢？

对于一个单位圆，在其内部做一个正 nnn 边形，要求每个顶点都在圆上，并且有一个顶点在 (1,0)(1,0)(1,0) 位置（因为有一个单位根一定是 111），那么这个 nnn 边形的每一个顶点就是一个单位根。

比如说当 n=6n=6n=6 时，图长这个样子：

其中，wnw_nwn 表示 nnn 次单位根，上面标记了 w6w_6w6 的点，都是单位根。

我们不妨设 (0,0)→(1,0)(0,0)\to(1,0)(0,0)→(1,0) 这条线段逆时针旋转遇到的第一条线段所对应的点为 wnw_nwn 的一次方，也就是 wn1w_n^1wn1，剩下的点依次为 wn2,wn3...w_n^2,w_n^3...wn2,wn3... 也就是：

能这样表示，那么自然也满足幂的运算规则：xa×xb=xa+bx^a \times x^b=x^{a+b}xa×xb=xa+b，在这里也就是 wna×wnb=wna+bw_n^a \times w_n^b=w_n^{a+b}wna×wnb=wna+b。

证明如下：

我们设原点到 wn1w_n^1wn1 这条线段与 xxx 轴的夹角为 α\alphaα，发现根据上面的编号规则，其实随着上标的增加，这个夹角也在增加，那么 wna+bw_n^{a+b}wna+b 的夹角大小就是 a×α+b×αa\times \alpha +b\times \alphaa×α+b×α。

根据上面的坐标的另外一种表示法，可以得到 wna+bw_n^{a+b}wna+b 的坐标为：(cos⁡(a×α+b×α)×1,sin⁡(a×α+b×α)×1)(\cos(a\times \alpha +b\times \alpha) \times 1,\sin(a\times \alpha +b\times \alpha) \times 1)(cos(a×α+b×α)×1,sin(a×α+b×α)×1)，也就是 (cos⁡(a×α+b×α),sin⁡(a×α+b×α))(\cos(a\times \alpha +b\times \alpha),\sin(a\times \alpha +b\times \alpha))(cos(a×α+b×α),sin(a×α+b×α))。

考虑两角和公式，有：
cos⁡(a×α+b×α)=cos⁡(a×α)cos⁡(b×α)−sin⁡(a×α)sin⁡(b×α)sin⁡(a×α+b×α)=sin⁡(a×α)cos⁡(b×α)+cos⁡(a×α)sin⁡(b×α)\cos(a\times \alpha + b\times \alpha)=\cos(a\times \alpha)\cos(b\times \alpha)-\sin(a\times \alpha)\sin(b\times \alpha)\\ \sin(a\times \alpha + b\times \alpha)=\sin(a\times \alpha)\cos(b\times \alpha)+\cos(a\times \alpha)\sin(b\times \alpha) cos(a×α+b×α)=cos(a×α)cos(b×α)−sin(a×α)sin(b×α)sin(a×α+b×α)=sin(a×α)cos(b×α)+cos(a×α)sin(b×α)

这便是 wna+bw_n^{a+b}wna+b 的坐标。

那么再看 wnaw_n^awna 乘 wnbw_n^bwnb 到底乘出了个什么东西。

先写出他们的坐标：
wna:(cos⁡(a×α),sin⁡(a×α))wnb:(cos⁡(b×α),sin⁡(b×α))w_n^a:(\cos(a\times \alpha),\sin(a\times \alpha))\\ w_n^b:(\cos(b\times \alpha),\sin(b\times \alpha)) wna:(cos(a×α),sin(a×α))wnb:(cos(b×α),sin(b×α))

要记得，这两个坐标相当于他们的实部和虚部。按照复数的运算规则，乘出来的数的坐标跟 wna+bw_n^{a+b}wna+b 一样，于是得证。

由此可以得出，单位根有五大性质：

wn0=wnn=1w_n^0=w_n^n=1wn0=wnn=1
wn0,wn1,...,wnn−1w_n^0,w_n^1,...,w_n^{n-1}wn0,wn1,...,wnn−1 互不相同
wnk=w2n2kw_n^k=w_{2n}^{2k}wnk=w2n2k
wnk=−wnk+n2(k<n2)w_n^k=-w_n^{k+\frac n 2}~~(k< \frac n 2)wnk=−wnk+2n (k<2n)
∑i=0n−1(wnj−k)i={0(j≠k)n(j=k)\sum_{i=0}^{n-1} (w_n^{j-k})^i=\begin{cases}0~~(j\neq k) \\n~~(j=k) \end{cases}∑i=0n−1(wnj−k)i={0 (j=k)n (j=k)

性质1证明： 过于明显就不证了。

性质2证明： 过于明显依然不证。

性质3证明：

感性理解：下标从 nnn 变成了 2n2n2n，意味着单位圆被分成了原来的两倍份，这样的话每一份的大小只有原来的 12\frac 1 221，那么因为我取的份数变成了原来的两倍，两倍乘 12\frac 1 221 等于 111，也就是和原来一样。

理性证明：带入定义式即可。
w2n2k=cos⁡(360°2n×2k)+sin⁡(360°2n×2k)i=cos⁡(360°n×k)+sin⁡(360°n×k)i=wnkw_{2n}^{2k}=\cos(\frac {~~360\degree}{2n} \times 2k)+\sin(\frac {~~360\degree}{2n} \times 2k)i~~~\\ =\cos(\frac {~~360\degree}{n} \times k)+\sin(\frac {~~360\degree}{n} \times k)i\\ =w_n^k~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ w2n2k=cos(2n 360°×2k)+sin(2n 360°×2k)i =cos(n 360°×k)+sin(n 360°×k)i=wnk

性质4证明：

感性理解：转半圈就变成原来的相反数了。

理性证明：带入定义式即可。
−wnk+n2=−wnk×wnn2=−wnk×(cos⁡(360°n×n2)+sin⁡(360°n×n2)i)=−wnk×(cos⁡(180°)+sin⁡(180°)i)=−wnk×(−1+0)=wnk-w_n^{k+\frac n 2}=-w_n^k \times w_n^{\frac n 2}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ =-w_n^k\times (\cos(\frac {~~360\degree} {n} \times \frac n 2) +\sin(\frac {~~360\degree} n \times \frac n 2)i)\\ =-w_n^k\times (\cos(180\degree) +\sin(180\degree)i)~~~~~~~~~~~~~~~~~~~~~~~\\ =-w_n^k\times (-1+0)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ =w_n^k~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ −wnk+2n=−wnk×wn2n =−wnk×(cos(n 360°×2n)+sin(n 360°×2n)i)=−wnk×(cos(180°)+sin(180°)i) =−wnk×(−1+0) =wnk

性质5证明： 难以感性理解，直接上证明。

当 j≠kj\neq kj=k 时，考虑等比数列求和。
∑i=0n−1(wnj−k)i=1−(wnj−k)n1−wnj−k=1−(wnn)j−k1−wnj−k=1−11−wnj−k=0\sum_{i=0}^{n-1} (w_n^{j-k})^i=\frac {1-(w_n^{j-k})^n} {1-w_n^{j-k}}~~~~~~~~~~~~~~~~~~~\\ =\frac {1-(w_n^n)^{j-k}} {1-w_n^{j-k}}\\ =\frac {1-1} {1-w_n^{j-k}}~~~~~\\ =0~~~~~~~~~~~~~~~~~~~~ i=0∑n−1(wnj−k)i=1−wnj−k1−(wnj−k)n =1−wnj−k1−(wnn)j−k=1−wnj−k1−1 =0

当 j=kj=kj=k 时，有：
∑i=0n−1(wnj−k)i=∑i=0n−1(wn0)i=∑i=0n−11=n\sum_{i=0}^{n-1} (w_n^{j-k})^i=\sum_{i=0}^{n-1} (w_n^0)^i=\sum_{i=0}^{n-1} 1=n i=0∑n−1(wnj−k)i=i=0∑n−1(wn0)i=i=0∑n−11=n

于是得证。

务必记得这些性质，当你看不懂下面的柿子时，就回来看看吧。

弧度制

在 C++C++C++ 里面，你是无法输入一个度数的，于是你需要弧度制。

具体来说，弧度制就是用一个弧来表示一个度数。

对于每一个度数 x°x\degreex°，他有一个对应的弧长，这段弧长的长度就是单位圆中做一个度数为 x°x\degreex° 的圆心角，这个圆心角所对弧的长度。

打个比方，比如说 60°60\degree60° 对应的弧长。

单位圆的半径是 111，那么周长就是 2π2\pi2π，那么这段弧的长度就是 2π6=13π\frac {2\pi} 6=\frac 1 3 \pi62π=31π。

在 C++C++C++ 中，13π\frac 1 3 \pi31π 就相当于 60°60\degree60°。

那么怎么在 C++C++C++ 中求 π\piπ 呢？

我们注意到，单位圆周长的一半就是 π\piπ，也就是说，π\piπ 对应 180°180\degree180°，因为 cos⁡(180°)=−1\cos(180\degree)=-1cos(180°)=−1，所以，我们可以用 acos(−1)acos(-1)acos(−1) 来算出 π\piπ。（注：acosacosacos 也就是 cos⁡\coscos 的逆运算，C++C++C++ 中是有这个函数的）

快速傅里叶变换

回到上面的话题，FFTFFTFFT 去找复数点，并不是指找复数平面直角坐标系上的点，而是去找以复数作为横坐标的点，因为复数平面直角坐标系上的点表示的终究只是一个数，而我们的最终目的是要找若干个真正的点。

要知道，复数既然是一个数，那么将这个数带入到函数式里面，也肯定是能得到一个值的，现在我们要做的，就是求出将 wn0,wn1,...,wnn−1w_n^0,w_n^1,...,w_n^{n-1}wn0,wn1,...,wnn−1 带入函数 CCC 中的值。

但是我们并不知道函数 CCC 的解析式呀，怎么求呢？

其实我们只需要求出将 xxx 带入函数 AAA 中的值 y1y_1y1 以及带入函数 BBB 中的值 y2y_2y2，然后乘起来就能得到将 xxx 带入函数 CCC 中的值，也就是 y1×y2y_1 \times y_2y1×y2。

现在问题变成了如何快速求出函数 AAA 和函数 BBB 在 wn0,wn1,...wnn−1w_n^0,w_n^1,...w_n^{n-1}wn0,wn1,...wnn−1 处的取值。

因为这两个问题是一样的，我们下面只考虑如何求出函数 AAA 在 wn0,wn1,...wnn−1w_n^0,w_n^1,...w_n^{n-1}wn0,wn1,...wnn−1 处的取值。

我们称，函数在 wn0,wn1,...wnn−1w_n^0,w_n^1,...w_n^{n-1}wn0,wn1,...wnn−1 处的取值为这个函数的傅里叶变换。

我们设函数 AAA 的解析式为 ∑i=0n−1aixi\sum_{i=0}^{n-1} a_ix^i∑i=0n−1aixi，将这 nnn 项根据 xxx 的上标的奇偶性分类，变成：
A(x)=(a0+a2x2+a4x4+...)+(a1x+a3x3+a5x5+...)A(x)=(a_0+a_2x^2+a_4x^4+...)+(a_1x+a_3x^3+a_5x^5+...) A(x)=(a0+a2x2+a4x4+...)+(a1x+a3x3+a5x5+...)

设 A1(x)=a0+a2x+a4x2+...,A2(x)=a1+a3x+a5x2+...A_1(x)=a_0+a_2x+a_4x^2+...~,~A_2(x)=a_1+a_3x+a_5x^2+...A1(x)=a0+a2x+a4x2+... , A2(x)=a1+a3x+a5x2+...，那么有：
A(x)=A1(x2)+xA2(x2)A(x)=A_1(x^2)+xA_2(x^2) A(x)=A1(x2)+xA2(x2)

我们考虑将 nnn 个单位根分成两部分来做，一部分是 wnk(k<n2)w_n^k~~(k<\frac n 2)wnk (k<2n)，另一部分是 wnk+n2(k<n2)w_n^{k+\frac n 2}~~(k<\frac n 2)wnk+2n (k<2n)。

先将 wnkw_n^kwnk 带入：
A(wnk)=A1((wnk)2)+wnkA2((wnk)2)=A1(wn2k)+wnkAx(wn2k)A(w_n^k)=A_1((w_n^k)^2)+w_n^k A_2((w_n^k)^2)~~~~\\ =A_1(w_{\frac n 2}^k) +w_n^kA_x(w_{\frac n 2}^k) A(wnk)=A1((wnk)2)+wnkA2((wnk)2) =A1(w2nk)+wnkAx(w2nk)

再将 wnk+n2w_n^{k+\frac n 2}wnk+2n 带入：
A(wnk+n2)=A1((wnk+n2)2)+wnk+n2Ax((wnk+n2)2)=A1((−wnk)2)−wnkAx((−wnk)2)=A1(wn2k)−wnkAx(wn2k)=A1(wn2k)−wnkAx(wn2k)A(w_n^{k+\frac n 2})=A_1((w_n^{k+\frac n 2})^2) +w_n^{k+\frac n 2}A_x((w_n^{k+\frac n 2})^2)~~~~~~~~~\\ =A_1((-w_n^k)^2) -w_n^kA_x((-w_n^k)^2)\\ =A_1(w_n^{2k}) -w_n^kA_x(w_n^{2k})~~~~~~~~~~~~~\\ =A_1(w_{\frac n 2}^k) -w_n^kA_x(w_{\frac n 2}^k)~~~~~~~~~~~~~~~ A(wnk+2n)=A1((wnk+2n)2)+wnk+2nAx((wnk+2n)2) =A1((−wnk)2)−wnkAx((−wnk)2)=A1(wn2k)−wnkAx(wn2k) =A1(w2nk)−wnkAx(w2nk)

诶？上下两个柿子只差一个符号！

这说明，当我们求出 A1A_1A1 和 A2A_2A2 时，就可以求出 AAA 了。

于是问题从求 AAA 在 wn0,wn1,...,wnn−1w_n^0,w_n^1,...,w_n^{n-1}wn0,wn1,...,wnn−1 处的取值变成了求 A1A_1A1 和 A2A_2A2 在 wn20,wn21,...,wn2n2−1w_{\frac n 2}^0,w_{\frac n 2}^1,...,w_{\frac n 2}^{\frac n 2-1}w2n0,w2n1,...,w2n2n−1 处的取值。问题的规模缩小了一半，这样就可以愉快地分治了！

于是代码如下：

void fft(comp *f,int len)
{if(len==1)return;//假如只有常数项，直接返回comp A1[len/2+1],A2[len/2+1];for(int i=0;i<(len>>1);i++)//按照奇偶性分类A1[i]=f[i*2],A2[i]=f[i*2+1];fft(A1,len>>1);fft(A2,len>>1);//递归下去comp wn(cos(2*pi/len),sin(2*pi/len)),w(1,0);//求出wn1(即wn)，w用来做累乘，求出wnkfor(int i=0;i<(len>>1);i++,w=w*wn)f[i]=A1[i]+A2[i]*w,f[i+(len>>1)]=A1[i]-A2[i]*w;//求解
}

做完了？不不不，滚动条告诉你并没有这么简单。

我们现在只是求出了这些点，我们还没得出解析式呢。

不对，点都还没求出来，还差这么一步：

fft(a,up);fft(b,up);
for(int i=0;i<up;i++)
c[i]=a[i]*b[i];

我们现在求出的 ccc 数组，相当于函数 CCC 的傅里叶变换，怎么通过这个傅里叶变换得到函数 CCC呢？

快速傅里叶逆变换

其实和快速傅里叶变换差不多，快速傅里叶变换是函数在 wn0,wn1,...wnn−1w_n^0,w_n^1,...w_n^{n-1}wn0,wn1,...wnn−1 处的取值，而逆变换就是在 wn0,wn−1,...wn−n+1w_n^0,w_n^{-1},...w_n^{-n+1}wn0,wn−1,...wn−n+1 处的取值。

显然 wn−1=wnn−1w_n^{-1}=w_n^{n-1}wn−1=wnn−1，因为 wn−1=wn−1×1=wn−1×wnn=wnn−1w_n^{-1}=w_n^{-1} \times 1=w_n^{-1}\times w_n^n=w_n^{n-1}wn−1=wn−1×1=wn−1×wnn=wnn−1。

这个很好实现，只需要将上面代码里的 wn1w_n^1wn1 换成 wn−1w_n^{-1}wn−1 即可，相当于顺时针旋转。发现 wn1w_n^1wn1 和 wn−1w_n^{-1}wn−1 的横坐标相同，纵坐标互为相反数，所以只需要将上面的

wn(cos(2*pi/len),sin(2*pi/len))

改成

wn(cos(2*pi/len),-sin(2*pi/len))

回到正题

我们设 ddd 为 ccc（这里的 ccc 是上面求出的 ccc 数组，也就是函数 CCC 的傅里叶变换）的傅里叶逆变换，那么有：
dk=∑i=0n−1ci(wn−k)i=∑i=0n−1(∑j=0n−1Cj(wni)j)(wn−k)i=∑i=0n−1∑j=0n−1Cj(wnj)i(wn−k)i=∑i=0n−1∑j=0n−1Cj(wnj×wn−k)i=∑i=0n−1∑j=0n−1Cj(wnj−k)i=∑j=0n−1Cj∑i=0n−1(wnj−k)id_k=\sum_{i=0}^{n-1} c_i(w_n^{-k})^i~~~~~~~~~~~~~~~~~~~~~~~~~\\ =\sum_{i=0}^{n-1} (\sum_{j=0}^{n-1}C_j (w_n^i)^j)(w_n^{-k})^i\\ =\sum_{i=0}^{n-1} \sum_{j=0}^{n-1}C_j (w_n^j)^i(w_n^{-k})^i~~~\\ =\sum_{i=0}^{n-1} \sum_{j=0}^{n-1}C_j (w_n^j\times w_n^{-k})^i~~\\ =\sum_{i=0}^{n-1} \sum_{j=0}^{n-1}C_j (w_n^{j-k})^i~~~~~~~~~~\\ =\sum_{j=0}^{n-1}C_j \sum_{i=0}^{n-1} (w_n^{j-k})^i~~~~~~~~~~~\\ dk=i=0∑n−1ci(wn−k)i =i=0∑n−1(j=0∑n−1Cj(wni)j)(wn−k)i=i=0∑n−1j=0∑n−1Cj(wnj)i(wn−k)i =i=0∑n−1j=0∑n−1Cj(wnj×wn−k)i =i=0∑n−1j=0∑n−1Cj(wnj−k)i =j=0∑n−1Cji=0∑n−1(wnj−k)i

这时候我们需要掏出我们封尘已久的性质5！

性质5：∑i=0n−1(wnj−k)i={0(j≠k)n(j=k)\sum_{i=0}^{n-1} (w_n^{j-k})^i=\begin{cases}0~~(j\neq k) \\n~~(j=k) \end{cases}∑i=0n−1(wnj−k)i={0 (j=k)n (j=k)

带入得：
=∑j=0n−1Cj×{0(j≠k)n(j=k)=Ck×n=\sum_{j=0}^{n-1}C_j \times \begin{cases}0~~(j\neq k) \\n~~(j=k) \end{cases}\\ =C_k \times n =j=0∑n−1Cj×{0 (j=k)n (j=k)=Ck×n

诶嘿？

也就是说，我们求出 ccc 数组的快速傅里叶逆变换后，将每一位除以 nnn，就可以得到 CCC 函数了！

你难道又以为完了？不不不，还有个细节问题。

发现我们每次需要将当前的序列进行奇偶性分类，分成两份，然后再合并得到当前序列的答案，那么这就要求，当前序列的长度是偶数。进而要求，原序列的长度需要是 222 的某个幂。

所以我们需要找到第一个大于序列长度的 222 的幂，让它来做序列长度，比原来多的位置不用管，放 000 即可。

代码是这样的：

while(up<=n+m)up<<=1;//n是函数A的长度，m是函数B的长度
//因为我们需要让up大于函数C的长度，所以需要让up>n+m

然后先给出逆变换的代码：

void fft(comp *f,int len,int type)//type表示做的是变换还是逆变换
{if(len==1)return;comp A1[len/2+1],A2[len/2+1];for(int i=0;i<(len>>1);i++)A1[i]=f[i*2],A2[i]=f[i*2+1];fft(A1,len>>1,type);fft(A2,len>>1,type);comp wn(cos(2*pi/len),type*sin(2*pi/len)),w(1,0);//注意这里type要乘上纵坐标for(int i=0;i<(len>>1);i++,w=w*wn)f[i]=A1[i]+A2[i]*w,f[i+(len>>1)]=A1[i]-A2[i]*w;
}

模板题完整代码：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 3000010int n,m;
struct comp{double x,y;comp(double xx=0,double yy=0):x(xx),y(yy){}comp operator +(const comp b){return (comp){x+b.x,y+b.y};}comp operator -(const comp b){return (comp){x-b.x,y-b.y};}comp operator *(const comp b){return (comp){x*b.x-y*b.y,x*b.y+y*b.x};}
};
comp a[maxn],b[maxn],c[maxn];
int up=1;
const double pi=acos(-1);
void fft(comp *f,int len,int type)
{if(len==1)return;comp A1[len/2+1],A2[len/2+1];for(int i=0;i<(len>>1);i++)A1[i]=f[i*2],A2[i]=f[i*2+1];fft(A1,len>>1,type);fft(A2,len>>1,type);comp wn(cos(2*pi/len),type*sin(2*pi/len)),w(1,0);for(int i=0;i<(len>>1);i++,w=w*wn)f[i]=A1[i]+A2[i]*w,f[i+(len>>1)]=A1[i]-A2[i]*w;
}int main()
{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);while(up<=n+m)up<<=1;fft(a,up,1);fft(b,up,1);for(int i=0;i<up;i++)c[i]=a[i]*b[i];fft(c,up,-1);for(int i=0;i<=n+m;i++)printf("%d ",(int)(c[i].x/(double)up+0.5));
}

很抱歉的告诉你，FFTFFTFFT 还没完。

这份优秀的代码提交之后，你可以获得 808080 分的好成绩。

为啥？

发现在分治的每一层我们都开了两个新的 A1A_1A1 和 A2A_2A2 数组，也就是说，我们有 nlognnlognnlogn 的空间是开在栈里面的，这就直接导致了栈空间爆炸。

于是，当数据大起来时，除了 ACACAC 之外，你还可以获得 TLETLETLE 或 RERERE。

于是，有了下一节。

迭代FFT

我们不能每次都把它进行奇偶性分类，那么不如一次分到底，然后从下往上合并？

假设我们已经分完了，那么代码就是这样的了：

void fft(comp *f,int len,int type)
{for(int mid=1;mid<len;mid<<=1)//枚举一块的大小{comp wn(cos(pi/mid),type*sin(pi/mid));//求出单位根for(int block=mid<<1,j=0;j<len;j+=block)//block是两块的大小，因为我们需要将相邻两块合成一整块，j枚举所有整块{comp w(1,0);for(int i=j;i<j+mid;i++,w=w*wn)//i枚举这一整块中的左边那一块，i+mid也就是在枚举右边那一块{comp x=f[i],y=f[i+mid]*w;//合并两块f[i]=x+y;f[i+mid]=x-y;}}}
}

是不是更简单啦~

剩下只有一个问题，如何求出每个数最后会被分到的位置？

找规律大法好！

观察一下每个数和他最后的位置的二进制数的关系。（第一行是这个数的二进制数，第二行是这个位置的二进制数）

相信大家肯定动动脚指头就能发现，每个数的二进制数的翻转就是他最后的位置。

那么暴力求即可。

我猜你又以为完了。

事实上并不是。

这个东西还有 O(n)O(n)O(n) 的求法~

对于一个数 iii，我们将他的最后一位拿走，然后进行翻转，再把拿走的那一位放在第一位，就可以得到 iii 的翻转。这就是 O(n)O(n)O(n) 的求法，因为把他的最后一位拿走之后的翻转我们已经知道了，所以可以直接得到答案。

打个比方：对于一个数 181181181，他的二进制数是 101101011011010110110101。

那么先拿走最后一位，变成：
10110101011010 1011010

翻转，变成
01011010101101 0101101

将一开始拿走的 111 放到第一位变成：
1010110110101101 10101101

于是我们得到了 iii 的翻转。

怎么知道 iii 被拿走最后一位的那个数的翻转呢？

iii 被拿走最后一位，相当于 i2\frac i 22i，那么这个翻转就是 i2\frac i 22i 的翻转除以 222。

为什么要除以 222？

因为 i2\frac i 22i 在翻转的时候，最高位被补了 000，翻转之后这个 000 就跑到了最低位，所以我们要除以 222 把它去掉。

依然用上面的那个栗子：

拿走最后一位之后，iii 变成
10110101011010 1011010

最高位补 000，变成：
0101101001011010 01011010

翻转，变成：
0101101001011010 01011010

除以 222，变成：
01011010101101 0101101

这才是我们想要的。

那么问题来了，为什么 i2\frac i 22i 的最高位会被补 000 ？

我们发现，在那个图片里，001001001 翻转变成 100100100，而不是 111 翻转变成 111，这说明，所有数被翻转的位数是固定的，不是根据这个数本身的位数而定的。而这个固定的位数，就是最大的那个数的位数。

于是我们可以得到递推式：

for(int i=1;i<up;i++)//第0位做了跟不做没区别
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
//r[i>>1]>>1 也就是去掉最低位后的翻转
//(i&1) 得到最低位，左移(l-1)位之后就是最高位了，然后与(r[i>>1]>>1)合并

模板题 ACACAC 代码：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define maxn 4000010const double pi=acos(-1);
int n,m;
struct comp{double x,y;comp(double xx=0,double yy=0):x(xx),y(yy){}comp operator +(const comp b){return (comp){x+b.x,y+b.y};}comp operator -(const comp b){return (comp){x-b.x,y-b.y};}comp operator *(const comp b){return (comp){x*b.x-y*b.y,x*b.y+y*b.x};}
};
comp a[maxn],b[maxn];
int r[maxn];
void fft(comp *f,int len,int type)
{for(int i=1;i<len;i++)//这一步下面再解释if(i<r[i])swap(f[i],f[r[i]]);for(int mid=1;mid<len;mid<<=1){comp w(cos(pi/mid),type*sin(pi/mid));for(int block=mid<<1,j=0;j<len;j+=block){comp p(1,0);for(int i=j;i<j+mid;i++,p=p*w){comp x=f[i],y=f[i+mid]*p;f[i]=x+y;f[i+mid]=x-y;}}}
}
void FFT(comp *p1,comp *p2,int len)
{fft(p1,len,1);fft(p2,len,1);for(int i=0;i<len;i++)p1[i]=p1[i]*p2[i];fft(p1,len,-1);
}int main()
{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 up=1,l=0;while(up<=n+m)up<<=1,l++;for(int i=1;i<up;i++)r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));FFT(a,b,up);for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/up+0.5));
}

发现代码里面有一个好像很奇怪的地方：

if(i<r[i])swap(f[i],f[r[i]]);

为什么需要判断 i<r[i]i<r[i]i<r[i] 呢？

因为 r[i]r[i]r[i] 是 iii 的翻转，假如有 r[a]=br[a]=br[a]=b，那么肯定有 r[b]=ar[b]=ar[b]=a。

假如我不加那个判断，就变成了：

swap(f[i],f[r[i]]);

那么当你枚举到 aaa 时，你会将 a,ba,ba,b 交换，枚举到 bbb 时，你又会将 b,ab,ab,a 交换，那么不就相当于没有变化？

所以，不妨设 a<ba<ba<b，那么加上那个判断之后，只有在枚举到 aaa 的时候会进行交换，而枚举到 bbb 时就不会交换了。

恭喜你，学完 FFTFFTFFT 了！

一点题表

ZJOI 2014 力题解
AH2017/HNOI2017 礼物题解