快速傅里叶想让我变换～天才们的分治多项式～

我看起来像是没事一样么……但是逃避又能怎样，就算害怕，我也要战斗。——《辉夜大小姐想让我表白~天才们的恋爱头脑战》

这东东并不属于中国计算机学会划定的提高组知识点考察范围，所以为了 提高组的小朋友可以退出这里了（除非你想多掉点头发）。

快速傅里叶变换

快速傅里叶变换（ ），是快速计算序列的离散傅里叶变换（ ）或其逆变换的方法。傅里叶分析将信号从原始域（通常是时间或空间）转换到频域的表示或者逆过来转换。主要用于解决一些多项式问题（卷积）问题，很大部分的现代技术，包括无线通讯，定位系统与广阔的信号处理领域相关的东西都基于该算法。

前置芝士有矩阵、离散傅里叶变换、时域到频域转换、原根、复数、微积分、欧拉公式等等，但是我一个都不会（真的，不是开玩笑）。

所以我们先从多项式乘法开始讲起。

多项式

多项式乘法

假设我们现在有两个多项式：

$$f(x)=x^2+3x+2\ \ \ \ \ \ \ g(x)=2x^2+1$$ $$h(x)=f(x)\cdot g(x)$$

然后现在我们重复应用乘法分配律来相乘多项式：

$$\begin{aligned}{} h(x) &= (x^2+3x+2)(2x^2+1)\\\\ &= x^2(2x^2+1) + 3x(2x^2+1)+2(2x^2+1)\\\\ &= 2x^4+x^2+6x^3+3x+4x^2+2\\\\ &= 2x^4+6x^3+5x^2+3x+2 \end{aligned}$$

非常好用，但是很慢，然而这就是乘法的本质，卷积。

卷积

如果现在用暴力来做多项式乘法，很轻松地得到以下方法：

设：

$$h(x)=f(x)\cdot g(x)$$

那么可以表达成：

$$\begin{aligned} &f[x]=\sum_{i=0}^{n}f[i]x^i\\\\ &g[x]=\sum_{i=0}^{n}f[i]x^i\\\\ &h[x]=\sum_{i=0}^kf[i]\cdot g[k-i] \end{aligned}$$

也就是想要得到 ，就需要 的第 项和 的第 项。

同样也可以写成这样：

$$h[x]=\sum_{i+j=k}f[i]\cdot g[j]$$

多项式乘法就是加法卷积，然后因为后面时间复杂度会噼里啪啦地炸开来所以就不讲了（擦汗）。

系数表达

然后有一个中朋友和大朋友能想到但是小朋友不能想到的想法就是将这个多项式的系数对应到一个序列里去：

$$\begin{aligned}{} &f(x)=x^2+3x+2 &\rightarrow \ \ \ &f=[2,3,1]\\\\ &g(x)=2x^2+1 &\rightarrow \ \ \ &g=[1,0,2]\\\\ &h(x)=2x^4+6x^3+5x^2+3x+2 &\rightarrow \ \ \ &h=[2,3,5,6,2]\\\\ \end{aligned}$$

列表中第 个系数对应到多项式的第 次项，这就是多项式的系数表示。

如果存在：

$$\begin{aligned}{} & f(x)=a_0+a_1x+a_2x^2+\cdots+a_dx^d\\\\ & g(x)=b_0+b_1x+b_2x^2+\cdots+b_dx^d\\\\ & h(x)=c_0+c_1x+c_2x^2+\cdots+c_{2d}x^{2d} \end{aligned}$$

这样子 里的每项系数都会去乘一边 里的每一项系数，时间复杂度就是是糟糕的 。

我们会环着这片海走下去。或许有一天，我们会真正回到起点

点值表示

一个真正不错的想法是点值表示。

次数为 的任何多项式都可以唯一地由 个点表示。

如果想要证明，和该点集对应的有且只有一组系数，如果我们实际计算经过这些点的多项式（这就是点值表示）：

$$\lbrace (x_0,f(x_0)),(x_1,f(x_1)),(x_d,f(x_d))\rbrace$$ $$\begin{aligned}{} & f(x)=f_0+f_1x+f2x^2+···+f_dx^d\\\\ & f(x_0)=f_0+f_1x_0+f_2x_0^2+\cdots+f_dx^d_0\\\\ & f(x_1)=f_0+f_1x_1+f_2x_1^2+\cdots+f_dx^d_1\\\\ & \qquad\qquad\qquad\qquad\ \ \vdots\\\\ & f(x_d)=f_0+f_1x_d+f_2x_d^2+\cdots+f_dx^d_d\\\\ \end{aligned}$$

然后我们就得到了一组方程式，就可以将其写成矩阵向量乘积来帮助分析（具体这个东西我也不是很清楚）：

$$\begin{bmatrix} f(x_0)\\\\ f(x_1)\\\\ \vdots\\\\ f(x_d) \end{bmatrix} = \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^d\\\\ 1 & x_1 & x_1^2 & \cdots & x_1^d\\\\ \vdots & \vdots & \vdots & \ddots & \vdots\\\\ 1 & x_d & x_d^2 & \cdots & x_d^d \end{bmatrix} \begin{bmatrix} f(x_0)\\\\ f(x_1)\\\\ \vdots\\\\ f(x_d) \end{bmatrix}$$

这里就不多加赘述了。

而我们将多项式的系数表示转换为点值表示，这个炫酷的魔术，就是闻名中外的快速傅里叶变换。

多项式求值

设我们现在有一个多项式：

$$f(x)=f_0+f_1x+f_2x^2+\cdots+f_dx^d$$

选取 个点，且满足 ，最简单的方法就是随机选取 个 坐标，然后计算对应的 坐标，但是当我们解构该过程时，又要使用到原多项式来求出 坐标，这样时间复杂度是 的，因为 ，兜兜转转又回到了起点，所以我们要找到优化的方式。

假设现在取一函数为 ，令 ，那么我们在这条函数上找那些知道一个点的函数值就知道另一点的函数值的点。

显然因为偶函数的定义 ，我们得到几个正负对的点。

再扩展到 的情况上：

此时因为奇函数的性质 ，得到以上的点（奇函数和偶函数是必修一的东西，小朋友可以感性理解）。

那么我们利用此想法扩展到更一般的多项式，对于一个多项式：

$$f(x)=3x^5+2x^4+x^3+7x^2+5x+1$$

分解这个多项式：

$$\begin{aligned} &f(x)=(2x^4+7x^2+1)+x(3x^4+x^2+5)\\\\ &f_e(x^2)=(2x^4+7x^2+1)\\\\ &f_o(x^2)=x(3x^4+x^2+5)\\\\ &f(x)=f_e(x^2)+xf_o(x^2) \end{aligned}$$

这样可以减少我们正负点对的计算，而且还把点的数量减少到了 个：

$$\begin{aligned} &f(x_i)=f_e(x_i^2)+x_if_0(x_i^2)\\\\ &f(-x_i)=f_e(x_i^2)-x_if_0(x_i^2) \end{aligned}$$

然后我们就得到了一个时间复杂度为 的优秀算法，但是其中有一个问题，在递归的时候，因为每个都是平方值，所有点的值最终都为正，所以递归会中断：

$$\begin{aligned} &f_e(x^2):[f_0,f_2,\cdots,f_{n-2}]\ \ \ \ \ \ \ \ \ \ \ f_o(x^2):[f_1,f_3,\cdots,f_{n-1}]\\\\ &[x_1^2,x_2^2,\cdots,x_{\frac{n}{2}}^2] \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_1^2,x_2^2,\cdots,x_{\frac{n}{2}}^2] \end{aligned}$$

唯一可行的方法是扩展初始点的值域，所以我们引入一个不得不引入的绝对领域，复数。

复数

在了解复数之前，首先我们要提到一个对于小朋友和中朋友很降智的问题：

$$\sqrt{-1}=?$$

诶，等等，根号里面不是不可以是负数的吗？

对于很多大朋友及以上的朋友，答案就是它了：

$$\sqrt{-1}=i$$

所以：

$$i^2=-1$$

这玩意不是我定义的所以我也不知道（笑）。

我们把形如 的数称为复数，其中的 称为实部， 称为虚部，复数域是实数域的代数闭包，即任何复系数多项式在复数域中总有根，复数域也是目前已知最大的域。在复平面中， 表示实数， 轴（除原点外的数）表示虚数，从原点 到 的向量表示复数 ，模长就是原点到点 的距离，即 。所谓幅角，假设以逆时针为正方向，从 轴正半轴到已知向量的转角的有向角叫做幅角。因为复数可以被表示为向量所以也可以使用平行四边形法则，这一点可在下面的代码中看到：

这里给出复数的代码，这份是使用重载运算符的：

struct Complex { // 复数类型  
	Complex (double xx = 0, double yy = 0) {x = xx, y = yy};
	double x, y;
	Complex operator + (Complex const &B) const // 加法
	{
		return Complex (x + B.x, y + B.y);
	}
	Complex operator - (Complex const &B) const // 减法
	{
		return Complex (x - B.x, y + B.y);
	}
	Complex operator * (Complex const &B) const // 乘法
	{
		return Complex (x * B.x - y * B.y, x * B.y + y * B.x);
	}
	Complex operator / (Complex const &B) const // 除法
	{
		double t = B.x * B.x + B.y * B.y;
		return Complex ((x * B.x + y * B.y) / t, (y * B.x - x * B.y) / t);
	}
} a, b; 

当然也有模板类版的（由 友情支持，但是跑得很慢，所以不推荐使用）：

typedef complex <double> comp; // 重定义一下（要定义就定义，不定义拉倒） 

单位根

上面好像只是一个广告，现在我们步入正题，假设现在有一个三阶多项式：

$$f(x)=x^3+x^2-x-1$$

那么我们需要 个点来表示该多项式，而且这些点还必须是正负对，那么可以想象这么一张图：

为了递归的可行性，这两点也必须是正负对，所以 和 之间需要一个等价关系，然后在底层，将只有一个 的点：

那么如果我们现在用 来代替这些数：

为了让右边等于 ，所以右上的两个数答案就是 和 。

这组特殊的点就是四次单位根（诺，就最底下的 ）。

那么如果我们有一个 次多项式，我们将需要 ，但是选 为 的整次方更简单，所以我们选 ，那么很简单了，推广到 次多项式，选择 ，使得 为 的整次方，选择的点集就是 次单位根。

将 的解是 次单位根均分在一个单位圆上，间隔点之间的角度是 ，简洁地列出这些点的方法就是欧拉公式：

$$e^{ie}=\cos(\theta)+i\sin(\theta)$$

这里要提到的一个东西是单位圆，在这里の三角函数采用弧度制，一个整圆不再是 ，而是 ，那么我们将其分成 份。

而 次单位根的标准定义是：

$$\omega=e^{\frac{2\pi i}{n}}$$

如果听不懂了，没有关系，后面的你更听不懂，这样就简洁地定义了每一个 次单位根，比如（就是上面写的那些）：

$$\begin{aligned} &w^0=1\\\\ &w^1=e^{\frac{2\pi i}{n}}\\\\ &w^2=e^{\frac{2\pi(2)i}{n}}\\\\ &w^{n-1}=e^{\frac{2\pi(n-1)i}{n}} \end{aligned}$$

然后，当我们要计算 次单位根处求一个多项式的值时，我们实际要做的就是求 の对应函数值，然后是 ，以此类推，直到 ，对于 ，求出它们在单位圆（模长为 ）上的的横坐标和纵坐标，实际上只需要算出 乘 次就能得到 了， 的坐标为 。

至于欧拉公式，是这么个东西（反正我是看不懂，不重要这个不是我们的重点，请感性理解，但是接下来你必须使用到它）：

让我们回到最初的问题，为什么要在 次单位根处用递归算法求 ？

观察上图发现，第 个 的根 与 为一对，在递归中我们将这些点进行平方，然后将其代入偶数项和奇数项两个多项式：

$$f(x)\rightarrow[1,\omega^1,\omega^2,...,\omega^{n-1}]$$ $$f_e(x^2)\ and\ f_o(x^2)\rightarrow[1,\omega^2,\omega^4,...,\omega^{2(\frac{n}{2}-1)}]$$

给出用欧拉公式求单位根的代码：

Complex sav (cos (2 * Pi / n), sin (2 * Pi / n)), buf (1, 0);
for (int i = 0; i < n; i ++)
{
	w[i] = buf; // w 因为很像 ω 所以就可以被使用（好吧其实我也想用 ω 奈何不彳亍） 
	buf = buf * sav; // 这里要是简写你就翘辫子了
}

FFT

DFT

讲了那么多，让我们来理一理 的过程：

$$f(x):[f_0,f_1,...,f_{n-1}]$$ $$\omega=e^{\frac{2\pi i}{n}}:[\omega^0,\omega^1,...\omega^{n-1}]$$ $$n=1\Rightarrow f(1)$$ $$f_e(x^2):[f_0,f_2,...,f_{n-2}]\ \ \ \ \ \ \ \ \ f_o(x^2):[f_1,f_3,...,f_{n-1}]$$ $$[\omega^0,\omega^2,...,\omega^{n-2}]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [\omega^0,\omega^2,...,\omega^{n-2}]$$ $$y_e=[f_e(\omega^0),f_e(\omega^2),...,f_e(\omega^{n-2})]\ \ y_o=[f_o(\omega^0),f_o(\omega^2),...,f_o(\omega^{n-2})]$$

先处理当 的时候，我们要对一点求值，这也是递归的出口。

接下来，该多项式要调用 两遍，一次在只有偶数项的多项式，一次在只有奇数项的多项式，这样新的多项式就是原来多项式的一半，我们只需要用 次单位根来求这些多项式的值，将偶数多项式和奇数多项式对应的点值表示结合起来以获取 多项式的点值表示。

其中利用到了正负对的思想：

$$f(x_j)=f_e(x_j^2)+x_jf_o(x_j^2)$$ $$f(-x_j)=f_e(x_j^2)-x_jf_o(x_j^2)$$ $$j \in \lbrace0,1,...\frac{n}{2}-1\rbrace$$

然后因为单位根的性质，就会产生以下关系：

$$x_j=\omega^j$$ $$-\omega^j=\omega^{j+\frac{n}{2}}$$

如果你懵掉了的话，我们需要讲一个欧拉恒等式：

$$e^{i\pi}=-1$$

利用单位根的性质，我们就可以把前面的第二个等式修改成：

$$f(\omega^{j+\frac{n}{2}})=f_e(\omega^{2j}-\omega^jf_o(\omega^{2j}))$$

然后，在代码里处理的时候：

$$y_e[j]=f_e(\omega^{2j})\ \ \ \ \ y_o[j]=f_o(\omega^{2j})$$

最后一步了（坚持住啊，后面还有更难的在等着你！），我们要返回在 次单位根处求出的多项式 的值：

$$y=[f(\omega^0),f(\omega^1),...,f(\omega^{n-1})]$$

看，这不就结束了吗（真的结束了吗？），现在，我们放出简单实现的代码：

#include <cstdio>
#include <cmath>

using namespace std;

const double Pi = acos (-1); // 结论，可以提高精度，环境允许也可以使用 M_PI（要用就用不用拉倒） 
const int N = 1e6 + 1;
int n, m;

struct Complex { // 复数类型  
	Complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
	double x, y;
	Complex operator + (Complex const &B) const // 重载运算符 
	{
		return Complex (x + B.x, y + B.y);
	}
	Complex operator - (Complex const &B) const
	{
		return Complex (x - B.x, y - B.y);
	}
	Complex operator * (Complex const &B) const
	{
		return Complex (x * B.x - y * B.y, x * B.y + y * B.x);
	}
	Complex operator / (Complex const &B) const // 实际上其实是用不到的  
	{
		double t = B.x * B.x + B.y * B.y;
		return Complex ((x * B.x + y * B.y) / t, (y * B.x - x * B.y) / t);
	}
} f[N], sav[N];

void dft (Complex *f, int len) // 传入数组和长度  
{
	if (len == 1) // 边界 
		return ;
	Complex *fl = f, *fr = f + len / 2; // 注意这里对指针巧妙的食用，注意完了吗，好，我也不懂  
	for (int k = 0; k < len; k ++)
		sav[k] = f[k]; // 复制数组  
	for (int k = 0; k < len / 2; k ++)
	{
		fl[k] = sav[k << 1]; // 分成两半 
		fr[k] = sav[k << 1 | 1]; 
	}
	dft (fl, len / 2); // 向下递归 
	dft (fr, len / 2); // 同上  
	Complex tG (cos (2 * Pi / len), sin (2 * Pi / len)), buf (1, 0); // 每次都要算单位根因为长度不同  
	for (int k = 0; k < len / 2; k ++)
	{
		sav[k] = fl[k] + buf * fr[k]; // 这里详见下面的公式（如果你不记得我们前面讲了什么的话，等等我前面好像没讲，讲了吗好像没讲讲了吧好像我记得算了管它讲没讲反正重复几遍总没有问题的） 
		sav[k + len / 2] = fl[k] - buf * fr[k]; // 同上  
		buf = buf * tG; // 得到下一个单位根，个十千百万不要简写 
	}
	for (int k = 0; k < len; k ++)
		f[k] = sav[k]; // 把数组复制过来  
} 

int main ()
{
	scanf ("%d", &n);
	for (int i = 0; i < n; i ++) 
		scanf ("%lf", &f[i].x); // 初始虚部为 0 ，所以不用管 
	for (m = 1; m < n; m <<= 1); // 将其补到 2 的幂次形式，反正都是 0 ，所以小问题 
	dft (f, m); // 传入数组和长度 
	for (int i = 0; i < m; i ++)
		printf ("(%.4lf,%.4lf)\n", f[i].x, f[i].y); // 这个输出有点奇怪啊，为什么不是多项式格式的呢？ 
	return 0;
}

给出上面说的主要公式：

$$f(\omega_n^{k+\frac{n}{2}})=f_l(\omega_{\frac{n}{2}}^k)-\omega_n^kf_r(\omega^k_{\frac{n}{2}})$$ $$f(\omega^k_n)=f_l(\omega^k_{\frac{n}{2}})+\omega_n^k f_r(\omega^k_{\frac{n}{2}})$$

因为是分治的那种感觉，所以时间复杂度是还能算优秀的 。

IDFT

诶诶诶上面的话都说了，还没有结束呢，做人啊，不能拉屎拉一半，屁股让别人擦（好像有啥不对但是我不记得那句话怎么讲了，反正就是咱们还有一个烂摊子要收就对了）。

上面的东西求完以后，只是求出了一个答案的点值表达，而并非最终结果，最后一步，我们要将求出的点值表达式还原成多项式，最终完成乘法，这个过程叫做傅里叶逆变换（不妨叫它大小姐逆告白 www）

然后这里有一个非常有意思的结论（涉及单位根反演，也可以理解成范德蒙德矩阵求逆）：

把 中的 换成 ，做完以后除以 即可，所以其实可以把 改一改就能叫 了（记住这句话）。

现在假设多项式 在变换之后得到的点值数列为 ，即 。

所以有以下式子：

$$g[k]=\sum_{i=0}^{n-1}(w_n^k)^if[i]$$

结论即：

$$n\cdot f[k]=\sum^{n-1}_{i=0}(\omega^{-k}_n)^ig[i]$$

最终的一个转换是，其中 是点值表示， 是向量：

$$a_k=\frac{c_k}{n}$$

所以我们把代码改成这样：

void idft (Complex *f, int len) // 逆快速傅里叶变换 
{
	if (len == 1)
		return ;
	Complex *fl = f, *fr = f + len / 2; // 注意这里对指针巧妙的食用，注意完了吗，好，我差不多懂了 
	for (int k = 0; k < len; k ++)
		sav[k] = f[k]; // 转移数组  
	for (int k = 0; k < len / 2; k ++)
	{
		fl[k] = sav[k << 1]; // 奇偶打乱  
		fr[k] = sav[k << 1 | 1];
	}
	idft (fl, len / 2); // 递归下去 
	idft (fr, len / 2);
	Complex tG (cos (2 * Pi / len), -sin (2 * Pi / len)), buf (1, 0);
	                // 贺代码的时候 ↑别忘了 
	for (int k = 0; k < len / 2; k ++)
	{
		// 这里 buf 是 len 次反单位根的第 k 次 
		sav[k] = fl[k] + buf * fr[k]; // 见上面公式，可能有点远要爪巴一下子 
		sav[k + len / 2] = fl[k] - buf * fr[k]; // 不对啊，你前面难道没看懂吗？为什么还要翻上去呢  
		buf = buf * tG; // 得到下一个单位根，十万百万千万不要简写
	}
	for (int k = 0; k < len; k ++)
		f[k] = sav[k]; 
} 

好像没啥改变，只有一个负号。

但是因为常数太大了所以时间复杂度会膨的一下爆炸，应该是过不了模板题。

代码实现

最后这里给出最简单的代码实现，我们可以把 和 合起来变成一个我们今天的主题：

$$dft+idft\rightarrow fft$$

#include <cstdio>
#include <cmath>

using namespace std;

const double Pi = acos (-1); // 结论，可以提高精度，环境允许也可以使用 M_PI（要用就用不用拉倒） 
const int N = 1e6 + 5e5 + 1; // 洛谷上的数据有点东西，如果只开 1e6 会 re 两个点 
int n, m;

struct Complex { // 复数类型
	Complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
	double x, y;
	Complex operator + (Complex const &B) const // 重载运算符 
	{
		return Complex (x + B.x, y + B.y);
	}
	Complex operator - (Complex const &B) const
	{
		return Complex (x - B.x, y - B.y);
	}
	Complex operator * (Complex const &B) const
	{
		return Complex (x * B.x - y * B.y, x * B.y + y * B.x);
	}
	Complex operator / (Complex const &B) const // 实际上其实是用不到的  
	{
		double t = B.x * B.x + B.y * B.y;
		return Complex ((x * B.x + y * B.y) / t, (y * B.x - x * B.y) / t);
	}
} f[N << 1], p[N << 1], sav[N << 1];

void fft (Complex *f, int len, bool flag) // 传入数组和长度，和一个判断（逆或不逆） 
{
	if (len == 1) // 边界 
		return ;
	Complex *fl = f, *fr = f + len / 2; // 注意这里对指针巧妙的食用，注意完了吗，好，我完全懂了  
	for (int k = 0; k < len; k ++)
		sav[k] = f[k]; // 复制数组  
	for (int k = 0; k < len / 2; k ++)
	{
		fl[k] = sav[k << 1]; // 分成两半 
		fr[k] = sav[k << 1 | 1]; 
	}
	fft (fl, len / 2, flag); // 向下递归 
	fft (fr, len / 2, flag); // 同上
	Complex tG (cos (2 * Pi / len), sin (2 * Pi / len)), buf (1, 0); // 每次都要算单位根因为长度不同  
	if (!flag)
		tG.y *= -1;
	for (int k = 0; k < len / 2; k ++)
	{
		sav[k] = fl[k] + buf * fr[k]; // 详见……自己上去看吧我不想再说一遍了 
		sav[k + len / 2] = fl[k] - buf * fr[k]; // 什么？？这里你也要我再说一遍？？ 
		buf = buf * tG; // 得到下一个单位根  
	}
	for (int k = 0; k < len; k ++)
		f[k] = sav[k]; // 把数组复制过来 
} 

int main ()
{
	scanf ("%d %d", &n, &m);
	for (int i = 0; i <= n; i ++)  
		scanf ("%lf", &f[i].x); // 初始虚部为 0 ，所以不用管 
	for (int i = 0; i <= m; i ++)
		scanf ("%lf", &p[i].x);
	for (m += n ,n = 1; n <= m; n <<= 1); // 将其补到 2 的幂次形式，反正都是 0 ，所以小问题 
	fft (f, n, 1); // 传入数组和长度，和判断
	fft (p, n, 1);
	for (int i = 0; i < n; i ++)
		f[i] = f[i] * p[i];
	fft (f, n, 0);
	for (int i = 0; i <= m; i ++)
		printf ("%d ", (int) (f[i].x / n + 0.49));
	return 0;
}

提交记录：

所以我们要优化！

优化

减少计算次数

来看这一段代码：

sav[k] = fl[k] + buf * fr[k]; // 都看到优化了你不会还不懂吧 
sav[k + len / 2] = fl[k] - buf * fr[k]; // emmm…… 

复数做乘法的代价是很高的，但是我们却偏偏把 算了两遍！

所以可以改成这样：

Complex temp = buf * fr[k]; // 先算出来 
sav[k] = fl[k] + temp; // 大哥……都到拆开来的时候了……这时候你要是再不懂 
sav[k + len / 2] = fl[k] - temp; // 我是真滴没有办法了 

蝴蝶变换

在代码中我们使用了大量的数组拷贝，现在要通过一些歪门邪道来避免数组拷贝（标题看着很阴间但是其实还好啦）。

让我们用一张图来看看分治究竟在搞什么飞机：

$$0\ \ \ \ \ \ 1\ \ \ \ \ \ 2\ \ \ \ \ \ 3\ \ \ \ \ \ 4\ \ \ \ \ \ 5\ \ \ \ \ \ 6\ \ \ \ \ \ 7$$ $$0\ \ \ \ \ \ 2\ \ \ \ \ \ 4\ \ \ \ \ \ 6\ \ \ \ \ \ 1\ \ \ \ \ \ 3\ \ \ \ \ \ 5\ \ \ \ \ \ 7$$ $$0\ \ \ \ \ \ 4\ \ \ \ \ \ 2\ \ \ \ \ \ 6\ \ \ \ \ \ 1\ \ \ \ \ \ 5\ \ \ \ \ \ 3\ \ \ \ \ \ 7$$

我们要先求出最后一层的数组状况，然后才能向上合并。

然后其实你可以观察一下（反正我是没看出来），有一个小小的规律：

$$6(110)_2\ \ \ \ \ 3(011)_2$$ $$4(100)_2\ \ \ \ \ 1(001)_2$$

是不是看出来了？看不出来我也要揭秘了，再不揭秘花儿都要谢了。

可以看到 在 原来的位置，而且它们的二进制数组是反过来的，这就是蝴蝶变换。

如果使用 的算法，就太糟糕了，我们要使用递推 的方法！代码如下：

for (int i = 0; i < limit; i ++)
	tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? limit >> 1 : 0); // tr[i] 是二进制反转

这个过程可以类比数位 （不会的可以溜去《动态规划目录》康康）

我们可以把一个二进制的反转拆成其他部分和最后一位两部分来看，其他部分的反转就是 tr[i >> 1] >> 1​ ，然后判一下最后一位，如果是 就加上 limit >> 1 。

也可以这样写：

for (int i = 0; i < limit; i ++)
{
	tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << bit);
	if (i < tr[i]) // 防止交换  
		swap (a[i], a[tr[i]]);
}

迭代

因为所有位置都确定了所以可以使用迭代的写法（而且这是最经典的写法，为什么不去学 呢）：

#include <bits/stdc++.h>

using namespace std;

const double Pi = acos (-1); // 结论，可以提高精度，环境允许也可以使用 M_PI（要用就用不用拉倒） 
const int N = 1e6 + 5e5 + 1; // 洛谷上的数据有点东西，如果只开 1e6 会 re 两个点 
int n, m;

struct Complex { // 复数类型
	Complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
	double x, y;
	Complex operator + (Complex const &B) const // 重载运算符 
	{
		return Complex (x + B.x, y + B.y);
	}
	Complex operator - (Complex const &B) const
	{
		return Complex (x - B.x, y - B.y);
	}
	Complex operator * (Complex const &B) const
	{
		return Complex (x * B.x - y * B.y, x * B.y + y * B.x);
	}
	Complex operator / (Complex const &B) const // 实际上其实是不得不得不得不用不到的  
	{
		double t = B.x * B.x + B.y * B.y;
		return Complex ((x * B.x + y * B.y) / t, (y * B.x - x * B.y) / t);
	}
} f[N << 1], p[N << 1];

int tr[N << 1];

void fft (Complex *f, bool flag) // 传入数组和长度，和一个判断（不得不逆或不逆） 
{
	for (int i = 0; i < n; i ++)
		if (i < tr[i])
			swap (f[i], f[tr[i]]); 
	for (int p = 2; p <= n; p <<= 1)
	{
		int len = p >> 1;
		Complex tG (cos (2 * Pi / p), sin (2 * Pi / p));
		if (!flag)
			tG.y *= -1;
		for (int k = 0; k < n; k += p) // 枚举起始点
		{
			Complex buf (1, 0);
			for (int l = k; l < k + len; l ++)
			{
				Complex tt = buf * f[len + l];
				f[len + l] = f[l] - tt;
				f[l] = f[l] + tt;
				buf = buf * tG; 
			}
		} 
	}
} 

int main ()
{
	scanf ("%d %d", &n, &m);
	for (int i = 0; i <= n; i ++)  
		scanf ("%lf", &f[i].x); // 初始虚部为 0 ，所以不用管 
	for (int i = 0; i <= m; i ++)
		scanf ("%lf", &p[i].x);
	for (m += n ,n = 1; n <= m; n <<= 1); // 将其补到 2 的幂次形式，反正都是 0 ，所以小问题 
	for (int i = 0; i < n; i ++)
		tr[i] = (tr[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0);
	fft (f, 1); // 传入数组和判断
	fft (p, 1);
	for (int i = 0; i < n; i ++)
		f[i] = f[i] * p[i];
	fft (f, 0);
	for (int i = 0; i <= m; i ++)
		printf ("%d ", (int) (f[i].x / n + 0.49));
	return 0;
}

如果要用建议直接用这个板子，下面是提交记录：

三次变两次

我刚开始看到这个我记得不是遇到鬼的时候要三步变两步来着说……，其实是减少调用函数的次数。

根据：

$$(a+bi)\times(c+di)=a(c+di)\times bi(c+di)=ac-bd+adi+bci$$

假设我们要求：

$$f(x)\times g(x)$$

设复多项式：

$$p(x)=f(x)+g(x)i$$

那么实部就是 ，虚部为 。

则：

$$p(x)^2=(f(x)+g(x)i)^2=f(x)^2-g(x)^2+2f(x)\cdot g(x)\cdot i$$

发现 的虚部就是 ，也就是说求出 之后，把它的虚部除以 即可。

那么主函数的那块只要这样写就可以了：

fft (f, 1); // 传入数组和判断
for (int i = 0; i < n; i ++)
	f[i] = f[i] * f[i];
fft (f, 0);
for (int i = 0; i <= m; i ++)
	printf ("%d ", (int) (f[i].y / n / 2 + 0.49));

然后快了很多：

注意，如果两个多项式的值域相差太大就会卡精度。

当然求多项式还有多种方法，像 （快速数论变换）， （快速沃尔什变换）， （某知名公司，也可以将其叫成快速谭炜谭变换）等等……