快速傅里叶变换

多项式全家桶

运算法则	算法	时间复杂度
多项式乘法	快速傅里叶变换	$\Theta (n{\log }_{2}n)$
多项式求逆	倍增+快速数论变换	$\Theta (n{\log }_{2}n)$
多项式对数函数	求导+积分	$\Theta (n{\log }_{2}n)$
多项式指数函数	泰勒展开+牛顿迭代	$\Theta (n{\log }_{2}n)$
分治FFT卷积	分治FFT/多项式求逆	$\Theta (n{\log }_2^{2}n)/\Theta (n{\log }_{2}n)$

---

作者的话

% 　　建议拥有部分高考知识、全部高考知识后再学习，本文将尽量照顾初中生。
　　本章为多项式全家桶基础中的基础——多项式乘法。

---

分类

缩写	全称	作用	时间复杂度
DFT	离散傅立叶变换	时频域转换	$O(n^2)$
FFT	快速傅立叶变换	时频域转换$($有精度误差$)$	$O(\text{大常数}+n{\log }_{2}n)$
NTT/FNTT	快速数论变换	模意义下的时频域转换	$O(\text{小常数}+n{\log }_{2}n)$
MTT	任意模数的NTT	任意模意义下的时频域转换	$O(n{\log }_{2}n)$
FWT	快速沃尔什变换	快速集合卷积	$O(\text{不定})$
FMT	快速莫比乌斯变换	逆莫比乌斯反演？	$O(\text{不定})$

---

快速傅里叶变换

前置技能

% 　　本文不包含但必不可少的前置技能：

1. sin函数，cos函数，弧度制
2. 函数
3. 笛卡尔坐标系(也叫做平面直角坐标系)
4. $\sum$ 的性质

卷积

% 　　什么是卷积？卷积是定义在函数上的运算。下面是百度百科的定义。

% 　　卷积是两个变量在某范围内相乘后求和的结果。如果卷积的变量是序列 $x(n)$ 和 $h(n)$，则卷积的结果为一个函数 $$y(n)=\sum _{i=-\infty }^{\infty }x(i)h(n-i)=x(n)\ast h(n)$$

% 　　出去的童鞋不要被百科拐跑了，百科哪有这里的文章没节操是不是…
　　修信号与系统的同学可能见过下面这段话：

% 　　给定义在 $(-\infty ,+\infty )$ 上的函数 $f_1(t)$ 与 $f_2(t)$ ，称由含参变量 $t$ 的广义积分所确定的函数$$g(t)={\int }_{-\infty }^{+\infty }{f}_1(\tau )f_2(t-\tau )d\tau$$　　为函数 $f_1(t)$ 与 $f_2(t)$ 的卷积，记为：$$g(t)=f_1(t)\ast f_2(t)$$

% 　　这里最简单的一句话概括：多项式乘法实际上是多项式系数向量的卷积。还是很晕？举个栗子吧。
　　设有两个多项式 $f=x^2+5x+4,g=3x^2-5x+7$。 则他们的卷积 $$h=f\times g=3x^4+10x^3-5x^2+15x+28$$　　其实是向量 $\vec{f}=(1,5,4),\vec{g}=(3,-5,7)$ 的卷积，为 $$\vec{h}=\vec{f}\otimes \vec{g}=(3,10,-5,15,28)$$　　根据多项式乘法的计算方法，可得：$$h_n=\sum _{i=0}^n{f}_i\cdot g_{n-i}$$　　现在再回去看百度百科的词条(就看我节选的那一段)，懂了吗？

系数表示法

% 　　设 $f(x)$ 为一个 $n$ 次函数，则其解析式为一个关于 $x$ 的 $n$ 次多项式。其解析式即为其系数表示法，也就是时域。
　　例：$f(x)=x^2+5x+4=(1,5,4)$
　　现在已知两个函数系数表示法，求这两个函数的卷积的系数表示法。显然地，需要 $O(n^2)$ 的时间。
　　设有两个函数 $f(x)=x^2+5x+4,g(x)=3x^2-5x+7$，则有
$$\begin{array}{rl}h(x)& =f(x)\otimes g(x)\\ & =(x^2+5x+4)\times (3x^2-5x+7)\\ & =3x^4+10x^3-5x^2+15x+28\end{array}$$

点值表示法

% 　　设 $f(x)$ 为一个 $n$ 次函数，因为 $n+1$ 个点确定一条 $n$ 次函数。因此可以选取函数上的 $n+1$ 个点来表示一个函数，称为点值表示法，也就是频域。例如：$$f(x)=\{(-1,0),(0,4),(1,10)\}$$　　如果已知两个函数在相同 $x$ 坐标上的点值表示法呢？求这两个函数的卷积的点值表示法。可以很负责任地说，只需要 $O(n)$ 的时间就可以完成。令
$$f(x)=\{(-1,0),(0,4),(1,10)\},g(x)=\{(-1,15),(0,7),(1,5)\}$$　　则有
$$\begin{array}{rl}h(x)& =f(x)\otimes g(x)\\ & =\{(-1,0),(0,4),(1,10)\}\otimes \{(-1,15),(0,7),(1,5)\}\\ & =\{(-1,0\times 15),(0,4\times 7),(1,10\times 5)\}\\ & =\{(-1,0),(0,28),(1,50)\}\end{array}$$　　可是不对啊，两个二次多项式的乘积应该是四次多项式啊，怎么能用三个点确定呢？ 没错，不能用三个点确定，那我们就用更多个点来确定。
$$\begin{array}{rl}h(x)& =f(x)\otimes g(x)\\ & =\{(-2,20),(-1,0),(0,4),(1,10),(2,18)\}\otimes \{(-2,29),(-1,15),(0,7),(1,5),(2,9)\}\\ & =\{(-2,20\times 29),(-1,0\times 15),(0,4\times 7),(1,10\times 5),(2,18\times 9)\}\\ & =\{(-2,580)(-1,0),(0,28),(1,50),(2,162)\}\end{array}$$ 　　但是一般情况下很少使用点值表示法，因此在使用这种方法之前需要进行转换，那转换的代价是多少？$$T(n)=\Theta (n)\times \Theta (n)=\Theta (n^2)$$　　好在J.W.库利和T.W.图基于1965年发明了快速傅里叶变换 (fast Fourier transform)。

向量

% 　　同时具有大小和方向的量。如下图中的向量OB记作 $\overrightarrow{OB}$

虚数和复数

% 　　定义虚数单位 ${\text{i}}^2=-1$(不用斜体以示区别)。则对于 $a,b\in \mathbb{R}$(实数域)$,b\ne 0$，形如 $z=a+b\text{i}$ 的数称为虚数。其中 $a$ 为实部，$b$ 为虚部。特别地，当且仅当 $a=0,b\ne 0$ 时，我们称 $z$ 为纯虚数。虚数和实数统称为复数。
　　在笛卡尔坐标系中，把 $x$ 轴当做实轴，把 $y$ 轴当做虚轴(单位长为 $\text{i}$)，这样的坐标系叫做复平面。
　　这上图中的向量对应了复数 $2+3\text{i}$ 。

struct complex{
	double x,y;//实部和虚部 
	complex (double xx=0,double yy=0) {
		x=xx,y=yy;
	}
};

复数的辐角　若复数 $z$ 所表示的向量与 $x$ 轴正半轴的夹角为 $\alpha$，则称 $\theta =\alpha +2k\pi ,k\in \mathbb{Z}$(整数域) 为复数 $z$ 的辐角，特别地，当 $\theta \in (-\pi ,\pi ]$ 时，我们称 $\theta$ 为复数 $z$ 的辐角主值，也称主辐角，记作 $\arg (z)$。如上图，$2+3i$ 的辐角为 $\theta +2k\pi ,k\in \mathbb{Z}$，辐角主值为 $\theta$。

复数的运算

% 　　加法运算：实部和虚部分别相加。例：
$$(3+5\text{i})+(8-4\text{i})=11+\text{i}$$　　减法运算：实部和虚部分别相减。例：
$$(3+5\text{i})-(8-4\text{i})=-5+9\text{i}$$　　乘法运算：多项式乘法。例：
$$\begin{array}{rl}(3+5\text{i})(8-4\text{i})& =3\times 8-3\times 4\text{i}+5\times 8\text{i}-5\text{i}\times 4\text{i}\\ & =24-12\text{i}+40\text{i}-20{\text{i}}^2\\ & =44+28\text{i}\end{array}$$

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

复数的模(长)　对于虚数 $z=a+b\text{i}$，定义 $z$ 的模(长)为 $$|z|=\sqrt{a^2+b^2}$$其几何意义是对应向量的模长，也是对应复平面上的点到原点的距离。

复数的三角形式　根据极坐标系的相关知识，在确定复数 $z$ 的辐角(主值) $\theta$ 与模长 $r=|z|$ 后复数有唯一解
$$z=r(\cos \theta +\text{i}\sin \theta )$$这便是复数的三角表示，显然有
$$a=r\cos \theta ,b=rsin\theta$$　　尽管复数的三角形式在加减法时显得十分笨拙，但复数的三角形式在乘除法方面上极具优势，一般地，对于两个复数$$x=r_1(\cos {\theta }_1+\text{i}\sin {\theta }_1),y=r_2(\cos {\theta }_2+\text{i}\sin {\theta }_2)$$根据三角恒等变换的知识，我们有$$\begin{array}{rl}x\times y& =r_1{r}_2(\cos {\theta }_1+\text{i}\sin {\theta }_1)(\cos {\theta }_2+\text{i}\sin {\theta }_2)\\ & =r_1{r}_2[(\cos {\theta }_1\cos {\theta }_2-\sin {\theta }_1\sin {\theta }_2)+\text{i}(\sin {\theta }_1\cos {\theta }_2+\cos {\theta }_1\sin {\theta }_2)]\\ & =r_1{r}_2[\cos ({\theta }_1+{\theta }_2)+\text{i}\sin ({\theta }_1+{\theta }_2)]\end{array}$$这说明两个复数相乘得到的复数模长等于两个复数模长之积，得到的复数辐角等于两个复数的辐角之和。特别地，复数 $z^n$ 的辐角为 $n\arg (z)$，模长为 $|z{|}^n$。
除法也同理，得到的复数主辐角相减，模长相除。
$$\begin{array}{rl}\frac{x}{y}& =\frac{r_1}{r_2}\frac{(\cos {\theta }_1+\text{i}\sin {\theta }_1)}{(\cos {\theta }_2+\text{i}\sin {\theta }_2)}\\ & =\frac{r_1}{r_2}\frac{(\cos {\theta }_1+\text{i}\sin {\theta }_1)(\cos {\theta }_2-\text{i}\sin {\theta }_2)}{(\cos {\theta }_2+\text{i}\sin {\theta }_2)(\cos {\theta }_2-\text{i}\sin {\theta }_2)}\\ & =\frac{r_1}{r_2}[\frac{(\cos {\theta }_1\cos {\theta }_2+\sin {\theta }_1\sin {\theta }_2)+\text{i}(\sin {\theta }_1\cos {\theta }_2-\cos {\theta }_1\sin {\theta }_2)}{{\cos }^2{\theta }_2+{\sin }^2{\theta }_2}]\\ & =\frac{r_1}{r_2}[\cos ({\theta }_1-{\theta }_2)+\text{i}\sin ({\theta }_1-{\theta }_2)]\end{array}$$

单位根(※)

代数(学)基本定理　任何复系数一元 $n$ 次多项式 方程在复数域上至少有一根($n⩾1$)。对这个定理的证明目前无法绕开高等数学，由这个定理，我们可以发现任何实系数一元 $n$ 次多项式方程在复数域上必有 $n$ 个复数根，只需根据代数基本定理分解因式分解出所有因式即可。

% 　　一般地，关于 $\omega$ 的方程 ${\omega }^n=1,n\in \mathbb{Z}$ 在实数域内有且仅有一个根 $1$，由代数基本定理，在复数域内，这个方程有 $n$ 个根，它们统称为 $n$ 次单位根，本文用 ${\omega }_n^i$ 来表示最小非负辐角第 $i+1$ 大的 $n$ 次单位根，可知 ${\omega }_n^0=1$。
　　事实上，这个方程在复数域内的$n$个根为 ${\omega }_n^i=\cos (i\cdot \frac{2\pi }{n})+\text{i}\sin (i\cdot \frac{2\pi }{n})$ ，因为其 $n$ 次乘方对应的辐角为 $n\cdot i\cdot \frac{2\pi }{n}=2i\pi$，故 $\arg \left(({\omega }_n^1{)}^n\right)=0$，$|({\omega }_n^1{)}^n|=1$，因而 $({\omega }_n^1{)}^n=1$，满足方程。
　　把它们画在复平面上可以发现，将单位圆 $n$ 等分，做 $n$ 个向量，那么这 $n$ 个向量对应复数即为 $n$ 次单位根，以 $1$ 为起点，逆时针分别为 ${\omega }_n^i$。
　　例如当 $n=8$ 时，在复平面上表示为(圆为单位圆)

% 　　规定这样的次序还可以带来一些便利，例如 ${\omega }_n^i=({\omega }_n^1{)}^i$，不妨再规定 ${\omega }_n^1={\omega }_n$，那么便可以使得 ${\omega }_n^i$ 既是第 $i$ 个 $n$ 次单位根，又表示 ${\omega }_n$ 的 $i$ 次方，同时将 ${\omega }_n$ 的上标由 $[0,n-1]$ 中的整数扩充到了 $\mathbb{Z}$。

单位根的性质(※)

性质1　$w_n^0=w_n^n=1$

性质2　$w_n^0,w_n^1,\cdots ,w_n^{n-1}$ 互不相同。

性质3　$w_{2n}^{2k}=w_n^k$
证　代入定义式即可，也可直接观察去掉单位圆上的奇数次单位根。
$$\begin{array}{rl}{w}_{2n}^{2k}=& \cos 2k\times \frac{2\pi }{2n}+\text{i}\sin 2k\times \frac{2\pi }{2n}\\ =& \cos k\times \frac{2\pi }{n}+\text{i}\sin k\times \frac{2\pi }{n}\\ =& w_n^k\end{array}$$性质4　$w_n^{k+\frac{n}{2}}=-w_n^k$
证　代入定义式即可，也可看出旋转半圈前后关于原点对称。
$$\begin{array}{rl}{w}_n^{k+\frac{n}{2}}=& w_n^k\times w_n^{\frac{n}{2}}\\ =& w_n^k\times \left[\cos \left(\frac{n}{2}\times \frac{2\pi }{n}\right)+\text{i}\sin \left(\frac{n}{2}\times \frac{2\pi }{n}\right)\right]\\ =& w_n^k\times (\cos \pi +\text{i}\sin \pi )\\ =& w_n^k\times (-1+0)\\ =& -w_n^k\end{array}$$性质5　$\begin{array}{r}\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=\{\begin{array}{l}0,k\ne j\\ n,k=j\end{array}\end{array}$
证　1. 当 $k\ne j$ 时，根据等比数列的求和公式，可得：
$$\begin{array}{rl}\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=& \frac{(w_n^{j-k}{)}^n-1}{w_n^{j-k}-1}=\frac{(w_n^n{)}^{j-k}-1}{w_n^{j-k}-1}=\frac{1-1}{w_n^{j-k}-1}=0\end{array}$$　　2.当 $k=j$ 时，可得：$\begin{array}{r}\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i=\sum _{i=0}^{n-1}1=n\end{array}$

---

正题

某些奇奇怪怪的读法。

1. $f\bar{a}f\bar{a}t\check{a}$
2. $FastFastTLE$

$\text{Cooley-Tukey}$ 算法

% 　　FFT 最常见的算法是 $\text{Cooley-Tukey}$ 算法，它的基本思路在 1965 年由 $\text{J.W.Cooley}$ 和 $\text{J.W.Tukey}$ 提出的，它是一个基于分治策略的算法。
% 　　显然地，一个 $n$ 次多项式可以被 $n+1$ 个点唯一确定。很恶心地，$\text{Cooley-Tukey}$ 算法在取点的时候不取整数点，不取有理数点，甚至还不取实数点！而是取复数点，$n$ 次单位根的 $0$ 到 $n-1$ 次幂。即使是这样，那时间复杂度仍然是 $\Theta (n^2)$ 的啊，而且复数的运算还比整数慢！别着急，来推推柿子。
　　设函数 $A(x)$ 的系数为 $(a_0,a_1,a_2,...,a_{n-1})$。则有：$$A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$$　　按照下标奇偶性分类，则有
$$A(x)=(a_0+a_2{x}^2+a_4{x}^4+...+a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+a_5{x}^5+...+a_{n-1}{x}^{n-1})$$
　　设$$\begin{array}{r}{A}_1(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{\frac{n}{2}-1}\\ A_2(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{\frac{n}{2}-1}\end{array}$$　　则有：$$A(x)=A_1(x^2)+x{A}_2(x^2)$$

% 　　这个时候，单位根终于派上用场了！
　　没错，我们把 $x=w_n^k(k<\frac{n}{2})$ 代入，得

$$\begin{array}{rl}A(w_n^k)& =A_1(w_n^{2k})+w_n^k{A}_2(w_n^{2k})\\ & =A_1(w_{\frac{n}{2}}^k)+w_n^k{A}_2(w_{\frac{n}{2}}^k)\end{array}$$

% 　　将 $x=w_n^{k+\frac{n}{2}}$ 代入，得：

$$\begin{array}{rl}A(w_n^{k+\frac{n}{2}})& =A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}{A}_2(w_n^{2k+n})\\ & =A_1(w_n^{2k}\ast w_n^n)-w_n^k{A}_2(w_n^{2k}\ast w_n^n)\\ & =A_1(w_n^{2k})-w_n^k{A}_2(w_n^{2k})\\ & =A_1(w_{\frac{n}{2}}^k)-w_n^k{A}_2(w_{\frac{n}{2}}^k)\end{array}$$　　发现了什么？只相差一个符号！有什么作用？那么当我们在枚举第一个式子的时候，我们可以 $O(1)$ 得到第二个式子的值。而且第一个式子的 $k$ 在取遍了 $[0,\frac{n}{2}-1]$ 的时候，第二个式子取遍了$[\frac{n}{2},n-1]$。因此原问题的规模缩小了一半，然后呢？
　　分治就好啦！下面是代码：

void fast_fast_tle(int limit,complex *a) {
	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];
	fast_fast_tle(limit>>1,a1);//分治
	fast_fast_tle(limit>>1,a2);
	complex Wn=complex(std::cos(2.0*Pi/limit),std::sin(2.0*Pi/limit)),w=complex(1,0);
	//Wn为单位根，w表示幂
	for(int i=0; i<(limit>>1); i++,w=w*Wn) //这里的w相当于公式中的k
		a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];//利用单位根的性质，O(1)得到另一部分
}

% 　　时间复杂度：$\Theta (\text{大常数}+nlo{g}_{2}n)$这个大常数是哪里来的呢？复数运算！复数的乘法是非常慢的。因此如果不是数据非常大，应尽量避免使用FFT。

点值乘法

% 　　经过了 $\text{Cooley-Tukey}$ 算法的折磨之后呢？然后直接乘就好了啊。

fast_fast_tle(n,a);
fast_fast_tle(n,b);
for(int i=0;i<=n;i++)
	a[i]=a[i]*b[i];

% 　　时间复杂度 $\Theta (n)$。完了？当然没有。还要逆回去。

快速傅里叶逆变换

% 　　设 $(y_0,y_1,...,y_{n-1})$ 是 $(a_0,a_1,a_2,\dots ,a_{n-1})$ 的傅里叶变换，即$(y_0,y_1,...,y_{n-1})$ 是 $(a_0,a_1,a_2,\dots ,a_{n-1})$ 在 $(w_n^0,w_n^1,\dots ,w_n^{n-1})$处的值。则有
$$\begin{array}{rl}{y}_i& =a_0+a_1{w}_n^1+a_2(w_n^1{)}^2+\cdots +a_n(w_n^1{)}^n\\ & =a_0+a_1{w}_n^1+a_2{w}_n^2+\cdots +a_n{w}_n^n\\ & =\sum _{j=0}^{n-1}{a}_j(w_n^i{)}^j\end{array}$$　　设 $(c_0,c_1,c_2,\dots ,c_{n-1})$ 是 $(y_0,y_1,...,y_{n-1})$ 在 $(w_n^0,w_n^{-1},\dots ,w_n^{-(n-1)})$ 处的取值。
% 　　求C序列同样用 $\text{Cooley-Tukey}$ 算法实现，只需要将只需要单位根变为 $w_n^{-1}=w_n^{n-1}$，相当于顺时针旋转。可以发现，在复平面上， $w_n^{n-1}$ 和 $w_n^1$ 的 $x$ 坐标相同， $y$ 坐标互为相反数，因此在代码中只需要把上面第9行的sin(2.0*Pi/limit)改成-sin(2.0*Pi/limit) 就可以了。
$$\begin{array}{rl}{c}_k& =\sum _{i=0}^{n-1}{y}_i(w_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j(w_n^i{)}^j)(w_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j(w_n^j{)}^i)(w_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j(w_n^j{)}^i(w_n^{-k}{)}^i)\\ & =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j(w_n^j{)}^i(w_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j(w_n^{j-k}{)}^i\\ & =\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i)\\ \therefore c_k& =\sum _{j=0}^{n-1}{a}_j(\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i)\end{array}$$　　根据性质5，有：

1. 当 $j\ne k$ 时，即 $j-k\ne 0$，此时 $\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i=0$，共有 $n-1$ 种情况
2. 当 $j=k$ 时，即 $j-k=0$，此时 $\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i=n$，共有 $1$ 种情况。
   $\therefore c_k=n\times a_k$

% 　　因而有
$$\frac{c_k}{n}=a_k$$　　然后呢？做完啦！先用 特别的 $\text{Cooley-Tukey}$ 算法求出 $(c_0,c_1,c_2,\dots ,c_{n-1})$，然后再除以 $n$ 就好啦！两种FFT代码有很多类似的部分，因此可以弄多一个参数，表示做的是那一种FFT。

void fast_fast_tle(int limit,complex *a,int type){
	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];
	fast_fast_tle(limit>>1,a1,type);
	fast_fast_tle(limit>>1,a2,type);
	//当type=1时，表示逆时针的FFT；当type=-1时，表示顺时针的FFT。 
	complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
	for(int i=0;i<(limit>>1);i++,w=w*Wn)
		a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}

% 　　执行完后还要除以 $n$，即：

for(int i=1;i<=limit;i++)
	a[i].x/=limit;

时空复杂度分析

% 　　令 $N$ 为大于 $n+m$ 的最小的 $2$ 的正整数次幂。
　　不难看出，FFT的时空复杂度均为 $O(\text{不小的常数}+Nlo{g}_{2}N)$
　　下面是完整的Luogu P3803代码

#include<cmath>
#include<cstdio>
const int MAXN=2*1e6+10;
const double Pi=std::acos(-1.0);
struct complex {
    double x,y;
    complex (double xx=0,double yy=0) {
        x=xx,y=yy;
    }
} a[MAXN],b[MAXN];
complex operator+(complex a,complex b) {
    return complex(a.x+b.x,a.y+b.y);
}
complex operator-(complex a,complex b) {
    return complex(a.x-b.x,a.y-b.y);
}
complex operator*(complex a,complex b) {
    return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
void fast_fast_tle(int limit,complex *a,int type){
    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];
    fast_fast_tle(limit>>1,a1,type);
    fast_fast_tle(limit>>1,a2,type);
    complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
    for(int i=0;i<(limit>>1);i++,w=w*Wn)
        a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}
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 limit=1;
    while(limit<=n+m) limit<<=1;
    fast_fast_tle(limit,a,1);
    fast_fast_tle(limit,b,1);
    for(int i=0; i<=limit; i++)
        a[i]=a[i]*b[i];
    fast_fast_tle(limit,a,-1);
    for(int i=0; i<=n+m; i++)
        printf("%d ",(int)floor(a[i].x/limit+0.5));
        //由于FFT涉及到浮点数运算，因此需要考虑精度误差。
    return 0;
}

% 　　于是乎，Luogu的提交结果无外乎有两种，一种是RE，另一种也是RE。
　　注意到空间复杂度中有$O(Nlo{g}_{2}N)$的空间是开在栈空间里的，因此可能会爆栈！所以有了下一小节。

迭代 $\text{Cooley-Tukey}$ 算法

% 　　观察一下原下标和新下标的二进制，发现了什么规律？
　　规律不是很好阐述，这里定义 $x$ 的 $l$ 位翻转为 $x$ 在长度为 ${\log }_{2}l-1$ 的二进制下的翻转¹。
　　举个例子：$1$ 的 $8$ 位翻转为 $4$。
　　$1$ 的长度为 ${\log }_{2}8-1=3$ 的二进制为 $001$。前后翻转为 $100$，十进制下即为 $4$。
　　可以发现，原下标为 $x$ 的新下标为 $x$ 的 $n$ 位翻转。设其为 $r[x]$。
　　若已经求出了 $r[i]$，则可以写出如下代码：

void fast_fast_tle(complex *A,int type) {
	for(int i=0;i<limit;i++)			//求出要迭代的序列
		if(i<r[i]) swap(A[i],A[r[i]]);
	for(int mid=1;mid<limit;mid<<=1) {	//待合并区间的长度(块的长度)
		complex Wn(std::cos(Pi/mid),type*std::sin(Pi/mid));//单位根
		for(int R=mid<<1,j=0;j<limit;j+=R) { //R是区间的右端点，j表示前已经到哪个位置了(枚举那一块)
			complex w(1,0);//幂
			for(int k=0;k<mid;k++,w=w*Wn) { //块内枚举
				complex x=A[j+k],y=w*A[j+mid+k];
				A[j+k]=x+y;
				A[j+mid+k]=x-y;
			}
		}
	}
}

% 　　那现在只剩下一个问题了：$r[i]$ 怎么求？

翻转序列

% 　　由于FFT的时间复杂度为 $O(Nlo{g}_{2}N)$。
% 　　因此如果求翻转序列的复杂度比 $O(Nlo{g}_{2}N)$ 还要高那就完了。

暴力

% 　　时间复杂度：$O(Nlo{g}_{2}N)$
% 　　空间复杂度：$O(n)$
% 　　代码复杂度：小
% 　　代码：下面的 $l={\log }_{2}limit$，实现上可以在求 $limit$ 的同时求出。

for(int i=1;i<=n-1;i++){
	int t=i,ret=0;
	for(int j=1;j<=l;j++){
		ret=(ret<<1)|(t&1);
		t>>=1;
	}
	r[i]=ret;
}

线性算法

% 　　想不到吧，这种东西还有线性算法。
% 　　下面以 $184$ 的 $256$ 位翻转为例子。
% 　　我们把一个数的二进制分成两部分，设这个数的二进制长度为 $l$ 位。
% 　　第一部分是前 $l-1$ 位，第二部分是最后一位。
$$1011100|0$$

% 　　$184$ 的 $256$ 位翻转等价于第二部分接上第一部分的 $128$ 位翻转。即
$$0|0011101$$

% 　　那第一部分的 $128$ 位翻转怎么求呢？等价于第一部分的 $256$ 位翻转右移 $1$ 位。
% 　　第一部分的 $256$ 位翻转
$$00111010$$

% 　　右移一位，得：
$$0011101$$

% 　　代码：

for(int i=1;i<limit;i++)
	r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));

时间复杂度：$O(n)$
空间复杂度：$O(n)$
代码复杂度：更小

最终代码

$\text{Accepted -O2}$ / 用时: $5111ms$ / 内存: $243752\text{KB}$

#include<bits/stdc++.h>
using namespace std;
const int maxn=2148576;
const double pi=std::acos(-1.0);
struct comp{
    double x,y;
    comp(double xx=0,double yy=0):x(xx),y(yy) {}
    friend comp operator+(const comp &x,const comp &y) {return comp(x.x+y.x,x.y+y.y);}
    friend comp operator-(const comp &x,const comp &y) {return comp(x.x-y.x,x.y-y.y);}
    friend comp operator*(const comp &a,const comp &b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
}a[maxn],b[maxn];
int limit=1,n,m,l=0,r[maxn];
void fft(comp *t,int ty){
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(t[i],t[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        comp wn(std::cos(pi/mid),ty*std::sin(pi/mid));
        for(int j=0,R=(mid<<1);j<limit;j+=R){
            comp w(1,0);
            for(int k=0;k<mid;k++,w=w*wn){
                comp x=t[j+k],y=w*t[j+k+mid];
                t[j+k]=x+y;
                t[j+k+mid]=x-y;
            }
        }
    }
}
int main(void)
{
    scanf("%d%d",&n,&m);
    while(limit<=n+m)
        limit<<=1,++l;
    for(int i=1;i<limit;i++)
        r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
    for(int i=0;i<=n;++i)
        scanf("%lf",&a[i].x);
    for(int i=0;i<=m;++i)
        scanf("%lf",&b[i].x);
    fft(a,1);
    fft(b,1);
    for(int i=0;i<limit;i++)
        a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=n+m;i++)
        printf("%d ",(int)(a[i].x/limit+0.5));
    return 0;
}

---

总结

% 　　恭喜完成一万字阅读。

---

题表(坑)

题目	来源	题解
A$\times$B Problem	Luogu 1919	构造
[MUTC2013]idiots	BZOJ3513	容斥原理

---

---

1. 其中 $l=2^x,x\in \mathbb{Z}$，$\log$的优先级比减法高。 ↩︎