概念引入

　　- 点值表示
　　　　对于一个$n - 1$次多项式$A(x)$，可以通过确定$n$个点与值（即$x$和$y$）来表示这唯一的$A(x)$

　　- 复数
　　　　对于一元二次方程
　　　　$$x^2 + 1 = 0$$
　　　　在实数范围内无解，那么我们将实数范围扩充，就得到了复数，再令$i$为该方程的解，即
　　　　$$i^2 = - 1$$
　　　　那么就定义$z = a + bi$的数为复数，则有
　　　　当$b = 0$时，$z$为实数
　　　　当$b \neq 0$时，$z$为虚数
　　　　当$a = 0 \&\& b \neq 0$时，$z$为纯虚数
　　　　其中，复数满足加法交换律、结合律、乘法分配率等

　　　　复数相乘：$z_1 = a_1 + b_1i, z_2 = a_2 + b_2i$，则$z_1 × z_2 = (a_1 + b_1i)(a_2 + b_2i) = (a_1a_2 - b_1b_2) + (a_1b_2 + b_1a_2)i$，它们相乘还是一个复数，在复平面上理解，就是模长相乘，幅角相加

　　　　共轭复数：当两个复数实部相同，虚部为相反数时，两个复数被称为共轭复数

　　　　对于一个复数$z = a + bi$，还可以在复数平面上用向量表示出来，即有$\overrightarrow{OZ}$一一对应每个$z = a + bi$，那么就有复数的模等于$\overline{Z}$，即$|z| = \sqrt{a^2 + b^2}$
　　　　复数直接比较大小没有意义，除非它是一个实数

$FFT$作用

　　那么，有了点值表示，对于多项式$A(x)$和$B(x)$相乘，就不需要$O(n^2)$，而只需要$O(n)$，因为$C(x_i) = A(x_i) * B(x_i)$，然后枚举$x_i$即可
　　于是现在的复杂度症结就在于将多项式转化成点值表示的$O(n^2)$了，于是就有$FFT$，可以实现在$O(n logn)$的时间内转化

离散傅里叶变换

　　于是傅里叶规定，点值中的$x$为模长为$1$的复数
　　至于为什么要取复数而不是实数，因为它有很神奇的性质
　　那么对于这$n$个复数的取法，取的是复平面上半径为$1$的一个圆，将其$n$等分后得到的点，令第$i$个点表示为$\omega_n^k(0 \le k \le n - 1)$，那么这个点在复平面中的表示即为$(cos\frac{k}{n}2\pi, sin\frac{k}{n}2\pi)$，那么根据复数乘法在复平面上的意义为模长相乘，幅长相加，即$\omega_n^k$相当于$(\omega_n^1)^k$，那么称$\omega_n^1$为单位根
　　我们就把这$\omega_n^0, \omega_n^1, …, \omega_n^{n - 1}$作为点值表示的$x$，称作离散傅里叶变换的结果
　　先给出关于傅里叶逆变换的结论：
　　　　- 将多项式$A(x)$的离散傅里叶变换结果作为系数代入多项式$B(x)$，再将每个离散傅里叶变换结果的倒数，即$\omega_n^0, \omega_n^{- 1}, …, \omega_n^{- (n - 1)}$作为$x$代入$B(x)$得到点值表示，那么这些表示除以$n$就得到了$A(x)$的原系数 [至于证明：不存在的]
　　这就是傅里叶变换神奇的性质

$FFT$

　　有了傅里叶变换，虽然多项式与点值表示相互转化已经很轻松了，但是复杂度仍然不理想，就有了快速傅里叶变换
　　结论：

  • $\omega_{2n}^{2k} = \omega_n^k$ [证明：代进原公式显然]
    • $\omega_n^{k + \frac{n}{2}} = - \omega_n^k$ [证明：关于复平面原点中心对称]
      对于多项式$A(x) = \sum\limits_{i = 0}^{n - 1} a_ix^i$，将它的奇偶项拆开，并将$x$转为$x^2$得到（此处先假定$n$为偶数）

$$
N(x) = a_0 + a_2x + a_4x^2 + … + a_{n - 2}x^{\frac{n}{2} - 1} \\
M(x) = a_1 + a_3x + a_5x^2 + … + a_{n - 1}x^{\frac{n}{2} - 1}
$$

　　 当然此时代入的是$N(x^2), M(x^2)$，则有

$$
A(x) = N(x^2) + xM(x^2)
$$
​ 　　此时就需要把$\omega_n^k$代入，分情况讨论：
　　 若$k < \frac{n}{2}$，则有
$$
\begin{aligned}
A(\omega_n^k) &= N(\omega_n^{2k}) + \omega_n^kM(\omega_n^{2k}) \\
&= N(\omega_{\frac{n}{2}}^k) + \omega_n^kM(\omega_{\frac{n}{2}}^k)
\end{aligned}
$$

​ 　　反之

$$
\begin{aligned}
A(\omega_n^k) &= N(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}}M(\omega_n^{2k + n}) \\
&= N(\omega_n^{2k} × \omega_n^n) + \omega_n^{k + \frac{n}{2}}M(\omega_n^{2k} × \omega_n^n)（这一步是通过复数乘法复平面上的意义化简的） \\
&= N(\omega_{\frac{n}{2}}^k) - \omega_n^kM(\omega_{\frac{n}{2}}^k)
\end{aligned}
$$

　　于是，我们只要求得$\{\omega_{\frac{n}{2}}^0, \omega_{\frac{n}{2}}^1, …, \omega_{\frac{n}{2}}^{\frac{n}{2} - 1}\}$，就可以得到$A(x)$的所有关于所有离散傅里叶变换结果的点值表示了，可用分治实现，复杂度$O(n logn)$

代码（分治FFT）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int MAXN = 2e06 + 10;

const double PI = acos (- 1.0);

int N, M;

struct mcomplex {
	double x, y;

	mcomplex () {}
	mcomplex (double fx, double fy) :
		x (fx), y (fy) {}

	mcomplex operator + (const mcomplex& p) const {
		return mcomplex (x + p.x, y + p.y);
	}
	mcomplex operator - (const mcomplex& p) const {
		return mcomplex (x - p.x, y - p.y);
	}
	mcomplex operator * (const mcomplex& p) const {
		return mcomplex (x * p.x - y * p.y, x * p.y + y * p.x);
	}
} ;
mcomplex omega (int n, int k, int inv) {
	return mcomplex (cos (2.0 * PI * (double) k / (double) n), 1.0 * inv * sin (2.0 * PI * (double) k / (double) n));
}
mcomplex A[MAXN], B[MAXN];

void FFT (mcomplex* a, int n, int inv) {
	if (n == 1)
		return ;
	mcomplex a1[n >> 1], a2[n >> 1];
	int m = n >> 1;
	for (int i = 0; i <= n; i += 2) {
		a1[i >> 1] = a[i];
		a2[i >> 1] = a[i + 1];
	}
	FFT (a1, m, inv);
	FFT (a2, m, inv);
	for (int i = 0; i < m; i ++) {
		mcomplex x = omega (n, i, inv);
		a[i] = a1[i] + x * a2[i];
		a[i + m] = a1[i] - x * a2[i];
	}
}

int getnum () {
	int num = 0;
	char ch = getchar ();

	while (! isdigit (ch))
		ch = getchar ();
	while (isdigit (ch))
		num = (num << 3) + (num << 1) + ch - '0', ch = getchar ();

	return num;
}

int main () {
	N = getnum (), M = getnum ();
	for (int i = 0; i <= N; i ++)
		A[i].x = (double) getnum ();
	for (int i = 0; i <= M; i ++)
		B[i].x = (double) getnum ();

	int n;
	for (n = 1; n <= N + M; n <<= 1);
	FFT (A, n, 1);
	FFT (B, n, 1);
	for (int i = 0; i <= n; i ++)
		A[i] = A[i] * B[i];
	FFT (A, n, - 1);
	for (int i = 0; i <= N + M; i ++) {
		if (i)
			putchar (' ');
		printf ("%d", (int) (A[i].x / n + 0.5));
	}
	puts ("");

	return 0;
}

蝴蝶操作

　　于是这样还是会超时，那么还需要优化
　　根据表格，有一个考眼力的性质
　　$|0 1 2 3|$
　　$|0 2 | 1 3|$
　　$|0 | 2 | 1 | 3|$
　　会发现每个数字的目标位置的二进制是原数的二进制翻转的结果，比如$1 = (01)_2, 2 = (10)_2$恰好是相对应的，于是就可以根据这个性质先将每个数排列到最终位置，再逐一合并

代码（迭代优化$FFT$）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int MAXN = (1 << 22);

const double PI = acos (- 1.0);

int N, M;

struct mcomplex {
    double x, y;

    mcomplex () {}
    mcomplex (double fx, double fy) :
        x (fx), y (fy) {}

    mcomplex operator + (const mcomplex& p) const {
        return mcomplex (x + p.x, y + p.y);
    }
    mcomplex operator - (const mcomplex& p) const {
        return mcomplex (x - p.x, y - p.y);
    }
    mcomplex operator * (const mcomplex& p) const {
        return mcomplex (x * p.x - y * p.y, x * p.y + y * p.x);
    }
} ;
mcomplex omega (int n, int k, int inv) {
    return mcomplex (cos (2.0 * PI * (double) k / (double) n), 1.0 * inv * sin (2.0 * PI * (double) k / (double) n));
}
mcomplex A[MAXN], B[MAXN];

int oppo[MAXN];
int limit;
void FFT (mcomplex* a, int inv) {
    for (int i = 0; i < limit; i ++)
        if (i < oppo[i])
            swap (a[i], a[oppo[i]]);
    for (int mid = 1; mid < limit; mid <<= 1) {
        mcomplex ome = mcomplex (cos (PI / (double) mid), inv * sin (PI / (double) mid));
        for (int n = mid << 1, j = 0; j < limit; j += n) {
            mcomplex x = mcomplex (1, 0);
            for (int k = 0; k < mid; k ++, x = x * ome) {
                mcomplex a1 = a[j + k], xa2 = x * a[j + k + mid];
                a[j + k] = a1 + xa2;
                a[j + k + mid] = a1 - xa2;
            }
        }
    }
}

int getnum () {
    int num = 0;
    char ch = getchar ();

    while (! isdigit (ch))
        ch = getchar ();
    while (isdigit (ch))
        num = (num << 3) + (num << 1) + ch - '0', ch = getchar ();

    return num;
}

int main () {
    N = getnum (), M = getnum ();
    for (int i = 0; i <= N; i ++)
        A[i].x = (double) getnum ();
    for (int i = 0; i <= M; i ++)
        B[i].x = (double) getnum ();

    int n, lim = 0;
    for (n = 1; n <= N + M; n <<= 1, lim ++);
    for (int i = 0; i <= n; i ++)
        oppo[i] = (oppo[i >> 1] >> 1) | ((i & 1) << (lim - 1));
    limit = n;
    FFT (A, 1);
    FFT (B, 1);
    for (int i = 0; i <= n; i ++)
        A[i] = A[i] * B[i];
    FFT (A, - 1);
    for (int i = 0; i <= N + M; i ++) {
        if (i)
            putchar (' ');
        printf ("%d", (int) (A[i].x / n + 0.5));
    }
    puts ("");

    return 0;
}

总结

　　这样的$FFT$可以在$O(n logn)$的时间内求出多项式乘法的各项系数，主要流程：将两个多项式分别转化成点值表示 $\dashrightarrow$ 通过点值表示将两个多项式合并 $\dashrightarrow$ 通过离散傅里叶逆变换将点值表示转化成系数表示，即得解

参考资料

[小学生都能看懂的FFT！！！]

[快速傅里叶变换(FFT)详解]

[复数——概念和代数运算]