Div1. 世界上最优雅的算法#

FFT起源

FFT的前身是DFT，可以简单看作是一堆OIer争先恐后对DFT算法进行优化的结果。美苏冷战期间，双方都对自己的核实力有所隐瞒，就等着某一天用自己的核导弹打对方个措手不及……在一次科学议会上，有人提出在苏联国境周边安装大量地震波传感器，将传感器的数据回收处理，分析是否出现了类似于核试验产生的震波从而判断苏联是否在进行秘密核试验。但是在当时用于频谱分析的DFT算法效率太过低下，传感器又必须在短时间内分析大量的频谱数据，因此逆境出人才——当时参与会议的其中一位科学家Tukey后来找到了程序员Cooley，后者在当时是一名操作ENIAC的程序员。二者分工明确，在1965年他们提出了FFT算法（当时叫做Cooley-Tukey算法），这个算法在后来被IEEE列入他们主编的“20世纪十大程序算法”之中。

FFT在信号处理领域受到了广泛的应用，直至今日它仍然被集成在大多数波谱分析仪器中作为底层算法之一。然而在信息学竞赛中，它被广大的OIer看中，拿去优化了高精度乘法的计算。它可以将朴素的高精度乘法的$\mathcal O(n^2)$复杂度大幅降低到$\mathcal O(n\log n)$，在处理大数据时，它将比Karatsuba乘法算法更为高效……

Div2. 系数表示法与点值表示法#

我们在初中时期学过一元二次方程的函数图像。一般情况下，我们都会用形如：$y=2x^2+x+1$的解析式形式来表示它。这种表示函数的方法就叫做系数表示法——自变量$x$的某次幂前乘一个系数得到的，一个$n$次多项式的系数表示法定义为：$y=\sum\limits_{i=0}^{n}a_ix^i$。

系数表示法可以直接通过$x$值求出$y$值，就好比生物的基因组排列，给出一个基因组，科学家就可以复原出整个生物的全貌，如果在各种函数绘制工具中输入解析式，它都唯一对应一个图像；但点值表示法就不是这样的了。它好比生物的某个性状，并且只是一个小得不能再小的特征点。比如给你两个特征：“白色毛发”和“部分呈现黑色”，不同人会做出不同的回答——动物学家会说“熊猫”、“斑马”、“美国短毛猫”一类的事物；阿宅会说：“伊雷娜”、“仆人”等等。加入其他限定词可能会使结果统一化，但是对于数学函数来说，你需要限定它的最高次数。如果你只给出了点$A(1,1)$，“求出经过点$A$的二次函数”这样的问题显然是不合适的，因为$n$次函数最少需要$n$个不同的点来唯一确定。

假如我们让一个三次函数$f(x)$和一个四次方程$g(x)$相乘求卷积，得到的函数$h(x)$将会是一个七次函数。如果我们需要将系数转换为点值，需要取定义域上七个不相同的$x$值进行计算，对于每次计算，要将$x$值以不同次幂代入至多七个算式中求出$y$值，显然时间复杂度是$\mathcal O(n^2)$的。那还不如不优化

为了减少这里的复杂度，我们借助一点高一所学的奇偶函数的内容~（没学过/没学懂没关系，奇函数相当于关于原点对称、偶函数相当于关于$y$轴对称，并且函数的取值范围也必须关于原点对称——假如$f(x)$能求出对应值，$f(-x)$也必须能求出对应值）。

假如该函数是一个偶次的幂函数，幂函数的定义是$f(x)=x^n,n\in\mathbb R$。如果此处$n$为偶数，整个函数就是一个偶函数、反之为奇函数。对于偶函数$e(x)$满足$e(x)=e(-x)$、奇函数$o(x)+o(-x)=0$。那么我们可以把原函数分解成多个偶函数和奇函数的和，例如：

函数$f(x)=6x^7+2x^5-x^4+x^3+2x^2+3$

经过一轮拆分：$f(x)=(6x^7+2x^5+x^3)+(-x^4+2x^2)+3$。根据奇偶函数的性质：奇函数+奇函数=奇函数、偶函数+偶函数=偶函数。因此左边括号整体组成的函数是奇函数、右边括号整体组成的函数是偶函数，常数单独拆出来：

这张图可以直观看出原函数（绿线）被拆分成偶项（红线）、奇项（蓝线）和常数项（橙线）后的样子。

对于奇项，代任意$x$，求出一个值后，第二个值直接取相反数；偶项则是$x$相反，结果相同。当然只拆一次是不够的，里面还有很多项呢，我们试着把原函数拆分成只有一项，以递归计算。我们对奇项和偶项继续进行拆分

令$f(x)=x(6x^6+2x^4+x^2)+(-x^4+2x^2)$，$f_o(x^2)=x(6x^6+2x^4+x^2),f_e(x^2)=-x^4+2x^2$。但是此时$f_e(x^2)$和$f_o(x^2)$的取值变成了平方，相对应地，继续拆分一次，括号里的$x^2$会变成$x^4$。因此取相反数$x$值不管用了（$x^2=(-x)^2=x\cdot x$），我们急切地想要找到就算平方后也是相反数的一对数据来求值。此时高阶思维娃找到了一个东西——复数域。

Div3. 单位根#

首先花一些篇幅来介绍复数的概念。我们学过：形如$x^2+1=0$这样的$\Delta<0$的一元二次方程是无解的。这并不完全正确——严格来说是在实数范围内无解，这个问题自从几千年前、方程发明之初就引来了无数数学家的疑问与探索，人们找来找去就是为了能找到一个阿拉伯数字，让它的平方等于$-1$。很可惜，找不到，因为传统意义上的阿拉伯数字体系建立在实数基础上，$\sqrt{-1}$又不在实数范围内，自然所有的尝试都会以失败告终。直到有一天，笛卡尔提出了虚数的概念，认为它是“不存在的数”、紧接着高斯使用虚数符号$i$来表示$\sqrt{-1}$的值，后面他又提出了“复平面”来将复数表示成平面上的一个向量。虚数这才有了立足之地。

而复数就像是实数和虚数的一次友好会面，复数域用字母$\mathbb C$表示，一个复数由实部和虚部组成。例如复数$3+9i$，它的实部是$3$、虚部是$9$。

紧接着，我们的FFT之旅就要来到下一个目的地——复平面。复平面类似于平面直角坐标系，只不过$x$轴变成了实数、$y$轴变成了虚数$i$。复平面的$x$轴因而称作实轴、$y$轴称作虚轴，刚才的复数$3+9i$可以看作复平面上的一个向量$(3,9)$。

比如说上一节遗留的式子，递归继续深入一层时，递归函数会变成$f_e(x^4)=\dots$。假如我们要让括号里的数是1，也就是这个函数在自变量为$1$时$x$的取值（事实上取成1恰好符合单位根的定义，于是就有简便的计算方式），我们就要解决$x^4=1$的问题，而且$x$还必须各不相同……两边同时开根号：$\pm x^2=\pm1$。展开成方程组就是：$\begin{cases}x^2=1&\dots1\\x^2=-1&\dots2\end{cases}$，根据$i$的数学定义，方程1解得$x_1=1,x_2=-1,x_3=-i^2$，方程2解得$x_4=i^2$。也就是：$\begin{cases}x_1=1\\x_2=-1\\x_3=-i^2\\x_4=i^2\end{cases}$。

每次递归，函数括号中的参数的幂就会乘2。幂在数值上就是$2^k,k\in\mathbb N_+$，解决$n$次多项式的拆分需要给出$n$个不同的$x$值，因而递归层数$p,p\in\mathbb N$满足如下关系：$2^p\geq n\rightarrow p\geq\log_2n$。因此原方程是七次方程，就需要递归三层，得到8个相异的$x$，从而计算出8个点！

引入单位根的概念：满足$z^n=1,n\in\mathbb N_+,z\in\mathbb C$的复数$z$叫做$n$次单位根。在复平面上，$n$次单位根将复平面上的单位圆$n$等分。回到上例中的四次单位根，它们在复平面上的分布如下图：

点$A$和点$C$恰好对应实数解$\pm1$，点$B$与点$D$对应复数解$\pm i^2$。

如果是八次（$2^3$）单位根，图像如下：

其中$E(\frac{\sqrt 2}{2},\frac{\sqrt 2}{2})$，对应复数$\omega_{n}^{1}=\frac{\sqrt 2}{2}+\frac{\sqrt 2}{2}i$，其余同理。由于$n$次单位根相当于将单位圆平分为$n$份，并且其中有一个解一定是$1$（上图$A$点）。联系到三角函数在单位圆中的表示，复平面单位圆上的点$(x,y)$一定满足$x=\cos\theta,y=i\sin\theta$。因此若$\omega_{n}^{0}=1$，则$\omega_{n}^{a}=\cos\frac{2a\pi}{n}+i\sin\frac{2a\pi}{n},a\in[0,n]$。

假如我们按逆时针方向，起始点$(1,0)$（实轴）从0开始编号，并将八次单位根和四次单位根的图像作比较，可以导出单位根的两条性质（$k\leq\frac{n}{2},k\in\mathbb N_+$）：

1. $\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$（对应点关于原点对称）
2. $\omega_{n}^{k}=\omega_{2n}^{2k}$（单位根次数和编号同乘2，坐标相同）

特殊地，$\omega_{n}^{0}=\omega_{n}^{n}=1$.

这下终于可以开始从一般例子出手推了。令$f(x)=a_0+a_1x+a_2x^2+a_3x^3+\dots+a_{2n-1}x^{2n-1}+a_{2n}x^{2n},n\in\mathbb N$。

按照奇偶项分为两组（这里可以直接通过$n$的奇偶性判断），折半次数，递归层数为1——令奇项函数为$f_o(x)$，偶项函数为$f_e(x)$。$f_o(x)=a_1+a_3x^2+a_5x^4+\dots+a_{2n-1}x^{2n-2}$（奇项提出一个$x$避免出现小数次幂）、$f_e(x)=a_0+a_2x^2+a_4x^4+\dots+a_{2n}x^{2n}$。将奇偶项的自变量设为$x^2$以表示出原式，即未折半的函数的未知数的次数：$f(x)=xf_o(x^2)+f_e(x^2)$。

函数的自变量来到了$x^2$，令$n=2$。该轮到单位根显神通了，令$k\leq\frac{n}{2},k\in\mathbb N_+$，此时的$x=\omega_{n}^{k}$，根据上述单位根的两条性质，推出如下式子：

$f(\omega_{n}^{k})=\omega_{n}^{k}f_o(\omega_{n}^{2k})+f_e(\omega_{n}^{2k})$，根据第二条性质，次数编号同除以2，得到：

$f(\omega_{n}^{k})=\omega_{n}^{k}f_o(\omega_{\frac{n}{2}}^{k})+f_e(\omega_{\frac{n}{2}}^{k})$；

如果此时$x=\omega_{n}^{k+\frac{n}{2}}$原式$=\omega_{n}^{k+\frac{n}{2}}f_o(\omega_{n}^{2k+n})+f_e(\omega_{n}^{2k+n})=\omega_{n}^{k+\frac{n}{2}}f_o(\omega_{n}^{2k}\cdot\omega_{n}^{n})+f_e(\omega_{n}^{2k}\cdot\omega_{n}^{n})$、根据第一条性质，$\omega_{n}^{k}=-\omega_{n}^{k+\frac{n}{2}}$以及特殊情况$\omega_{n}^{n}=1$，代入，根据第二条性质，次数编号同除以2，得到：

$f(\omega_{n}^{k+\frac{n}{2}})=-\omega_{n}^{k}f_o(\omega_{n}^{2k})+f_e(\omega_{n}^{2k})=-\omega_{n}^{k}f_o(\omega_{\frac{n}{2}}^{k})+f_e(\omega_{\frac{n}{2}}^{k})$。

正因此，当$\omega_{\frac{n}{2}}^{1},\omega_{\frac{n}{2}}^{2},\omega_{\frac{n}{2}}^{3}\dots\omega_{\frac{n}{2}}^{\frac{n}{2}-1},\omega_{\frac{n}{2}}^{\frac{n}{2}}$已知时，代入$f(x)$计算，那么当$x$取$\omega_{\frac{n}{2}}^{1+\frac{n}{2}},\omega_{\frac{n}{2}}^{\frac{n}{2}+2},\omega_{\frac{n}{2}}^{\frac{n}{2}+3}\dots\omega_{\frac{n}{2}}^{n-1},\omega_{\frac{n}{2}}^{n}$时的$f(x)$值也就明确了（根据性质第二条直接求出）。$f_o(x)$和$f_e(x)$就可以将问题规模缩小一半，递归的边界就是$n=1$。

Div4. IFFT与蝴蝶操作#

IFFT的中文名称是快速傅里叶逆变换。有些人比如刚开始接触FFT的我认为只要敲出FFT的代码，就可以解决任何函数卷积的问题了……错了，FFT主要是将系数表示法转换为点值表示法，相当于数学考试解析式求解题题干里给你的已知点坐标！卷积之路才刚刚过半，接下来介绍如何将相乘后得到的点值表示法重新转化为系数表示法。由于FFT是系数到点值的转换、这种方式是点值到系数的转换，因而称之快速傅里叶逆变换——IFFT。

我们的FFT算法已经将两个函数的卷积$h(x)=f(x)g(x)$按照点值形式求了出来，但是一般生活中大家还是习惯性使用解析式，也就是系数表示法表示一个函数，就需要借助IFFT将卷积的点值形式重新转化为系数形式。

假设一个函数的系数形式$v$在经历了一次FFT后变成了其点值形式$p$。根据FFT有：$p_k=\sum\limits_{i=0}^{n-1}v_i\cdot \omega_{n}^{ik}$。IFFT的作用就是将$p_k$变成$v_k$，下面是具体操作：

$v_k=\frac{1}{n}\sum\limits_{i=0}^{n-1}p_i\omega_{n}^{-ik}$

简记为：根倒求和、$n$倒乘积。

再加上，复数的乘积在复平面上有个奇妙的性质：复平面上一个向量可以看作是从实轴正半轴逆时针旋转特定的角度后得到的，这个角度称作向量的辐角，复数相乘时，得到的结果的辐角将是这两个向量的辐角之和、模长将是两个向量模长的乘积。简称辐角相加、模长相乘。除法作为乘法的逆运算，它在复平面上操作的结果是辐角相减、模长相除。因而单位根的倒数$\omega_{n}^{-k}$可以看作$1\div\omega_n^k$，也就是从实轴正半轴顺时针旋转$\omega_n^k$的辐角大小。计算$\omega$时就可以把原单位根关于实轴做一次对称变换，即虚部取相反数，也就是共轭。

对于代码的设计，尽管C++STL库提供了Complex<>复数类，但还是建议手写Complex。因为STL库的运行效率有时会很慢，会被卡掉，而且它的码量也不大，手写难度不高。强烈推荐手写结构体Complex。

递归Version：

1
#include <bits/stdc++.h>
2
#define N 6000010
3
using namespace std;
4

5
const double PI = acos(-1.0);
6

7
struct Complex {
8
  double real, imag;
9

10
  Complex(double r = 0.0, double i = 0.0) {
11
    real = r, imag = i;
12
  }
13
} f[N], g[N];
14

15
Complex operator +(const Complex &l, const Complex &r) {
16
  Complex res(l.real + r.real, l.imag + r.imag);
17
  return res;
18
}
19

20
Complex operator -(const Complex &l, const Complex &r) {
21
  Complex res(l.real - r.real, l.imag - r.imag);
22
  return res;
23
}
24

25
Complex operator *(const Complex &l, const Complex &r) {
26
  //(a+bi)*(c+di)=(ac-bd)+(ad+bc)i
27
  Complex res(l.real * r.real - l.imag * r.imag, l.real * r.imag + l.imag * r.real);
28
  return res;
29
}
30

31
Complex operator /(const Complex &l, const Complex &r) {
32
  //(a+bi)/(c+di)=(a+bi)(c-di)/(c^2+d^2)=[(ac+bd)+(bc-ad)i]/(c^2+d^2)
33
  Complex res((l.real * r.real + l.imag * r.imag) / (r.real * r.real + r.imag * r.imag), (l.imag * r.real - l.real * r.imag) / (r.real * r.real + r.imag * r.imag));
34
  return res;
35
}
36

37
void FT(int len, Complex *c, int type) {
38
  // type=1时为FFT，type=-1时单位根纵坐标取相反数，为IFFT
39
  if (len == 1) return;
40
  Complex c1[len >> 1], c2[len >> 1]; // c1为偶项系数，c2为奇项系数，数组大小除以2
41
  for (int i = 0; i <= len; i++) {
42
    // 按奇偶性分类存储
43
    if (i % 2 == 0) c1[i >> 1] = c[i]; // 紧凑存储，c1、c2数组下标除以2
44
    else c2[i >> 1] = c[i];
45
  }
46
  FT(len >> 1, c1, type); // 处理偶项
47
  FT(len >> 1, c2, type); // 处理奇项
48
  Complex omega = Complex(cos(2.0 * PI / len), type * sin(2.0 * PI / len)); // 计算单位根
49
  Complex k = Complex(1.0, 0.0);
50
  register Complex butterfly;
51
  for (int i = 0; i < (len >> 1); i++, k = k * omega) {
52
    butterfly = k * c2[i]; // 蝴蝶操作，记录下复数乘积（现场算太慢）便于调用
53
    c[i] = c1[i] + butterfly; // 单位根
54
    c[i + (len >> 1)] = c1[i] - butterfly; // 单位根性质，系数取反得到另一半
55
  }
56
}
57

58
void IFFT(int len, Complex *c) {
59
  FT(len, c, -1);
60
}
61

62
void FFT(int len, Complex *c) {
63
  FT(len, c, 1);
64
}
65

66
int main() {
67
  int n, m;
68
  cin>>n>>m;
69
  for (int i = 0; i <= n; i++) cin>>f[i].real;
70
  for (int i = 0; i <= m; i++) cin>>g[i].real;
71
  int len = 1;
72
  while (len <= n + m) len <<= 1; // 得到最小递归层数
73
  FFT(len, f);
74
  FFT(len, g);
75
  for (int i = 0; i <= len; i++) f[i] = f[i] * g[i]; // 点值相乘得到卷积点值
76
  IFFT(len, f); // 卷积点值转化为系数
77
  for (int i = 0; i <= n + m; i++) cout<<(int) (f[i].real / len + 0.5)<<' ';
78
  return 0;
79
}

时间复杂度：$\mathcal O(n\log n)$

这个代码足够通过模板题有些人说不行，只需把数组开大一点点（5e6起步）就可以了……但是这个代码其实还可以继续优化，但是这下就很棘手了，本来推来推去已经够烧脑了你还要让我优化？，考虑到递归算法有些低效，浪费了一些空间和执行效率。那么我们如何把它变成一个非递归版本呢？答案是迭代！

迭代FFT与二进制优化#

我们按照奇偶下标将原多项式分为了偶项和奇项，这样做的复杂度是$\mathcal O(n\log n)$。我们仔细观察一下抽出的系数之间下标的联系：

递归之前，系数下标：$0\;1\;2\;3\;4\;5\;6$

递归一层，系数下标：$0\;2\;4\;6\;|\;1\;3\;5$

递归两层，系数下标：$0\;4\;|\;2\;6\;\;\;|\;\;\;1\;5|\;3$

递归三层，系数下标：$0\;|4\;|\;2\;|\;6\;|\;1\;|\;5\;|\;3$

既然称作二进制优化，我们来看看第三层的数字写成二进制形式是怎样的：000 100 010 110 001 101 011

没有规律……别急，我们把数字倒过来念：000 001 010 011 100 101 110。诶？他们是单调递增的，再写成十进制就是：0 1 2 3 4 5 6。我们就找到规律了：递归结束时系数下标的二进制值等于原多项式的系数下标的颠倒二进制值

常规操作，我们使用进制转换来实现二进制逆序，这样的时间复杂度是$\mathcal O(\log n)$的：

1
int rev(int n) {
2
  int res = 0;
3
  int k = log2(n) + 1;
4
  for (int i = 1; i <= k; i++) {
5
    if (n & 1) res += pow(2, k - i);
6
    n >>= 1;
7
  }
8
  return res;
9
}

借鉴了这篇题解的逆序思路，逆序的时间复杂度是$\mathcal O(1)$的，处理数组的复杂度是$\mathcal O(n)$的：

1
int l,r[MAXN];
2
int limit=1;
3
void fast_fast_tle(complex *A,int type)
4
{
5
    for(int i=0;i<limit;i++)
6
        if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列
7
    for(int mid=1;mid<limit;mid<<=1)//待合并区间的中点
8
    {
9
        complex Wn( cos(Pi/mid) , type*sin(Pi/mid) ); //单位根
10
        for(int R=mid<<1,j=0;j<limit;j+=R)//R是区间的右端点，j表示前已经到哪个位置了
11
        {
12
            complex w(1,0);//幂
13
            for(int k=0;k<mid;k++,w=w*Wn)//枚举左半部分
14
            {
15
                 complex x=A[j+k],y=w*A[j+mid+k];//蝴蝶效应（应为“蝴蝶操作”）
16
                A[j+k]=x+y;
17
                A[j+mid+k]=x-y;
18
            }
19
        }
20
    }
21
}

以及r数组的读入：

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

当然，提高组不考，如果考到，基本上那个递归版本也够用了……