快速傅里叶变换
多项式
定义
不严谨定义：形如 \(f(x) = \sum \limits _{i=0}^{n} a_ix^i\) 的式子为多项式。
定义（from OIWiki）：对于求和式 \(\sum a_nx^n\)，如果是有限项相加，称为多项式，记作 \(f(x)=\sum \limits_{n=0}^m a_nx^n\)。
次数：对于多项式 \(F(x) = \sum \limits _{i=0}^{n} a_ix^i\)，该多项式的次数为 \(n\)。表示为 $\text{degree} (F) = n $
次数界：任何一个大于多项式次数的整数。
基本运算
系数表示法：
\[\mathbf{a} = [a_1, a_2,\cdots, a_n] ^\top \]加法
\[\mathbf{c} = \mathbf{a} + \mathbf{b} \\ c_i = \sum \limits_{i = 0}^{n} (a_i + b_i) \]时间复杂度 \(O(n)\)。
乘法
\[A(x) = \sum \limits_{i = 0}^n a_ix^i \ \ B(x) = \sum \limits_{i = 0}^n b_ix^i \\ \text{degree} (C) = \text{degree} (A) + \text{degree} (B) \\ c_i = \sum \limits_{j = 0}^i a_j b_{i - j} \]记为 \(\mathbf{c} = \mathbf{a} \otimes \mathbf{b}\)
时间复杂度 \(O(n^2)\)。
求值
秦九韶算法：
\[\begin{aligned} f(x) ={}& a_nx^n+a_{n-1}x^{n-1}+\cdots a_1x+a_0 \\ ={}& (a_nx^{n-1}+a_{n-1}x^{n-2}+\cdots a_1)x+a_0 \\ ={}& ((\cdots(a_nx+a_{n-1})x+a_{n-2})x+\cdots))x+a_0 \end{aligned} \]时间复杂度 \(O(n)\)。
点值表示法：
\[\{(x_0, A(x_0)), (x_1, A(x_1)), \cdots(x_n, A(x_n))\} \]加法
\[C(x_i) = A(x_i) + B(x_i) \]时间复杂度 \(O(n)\)。
乘法
\[C(x_i) = A(x_i)B(x_i) \]时间复杂度 \(O(n)\)。
求值
拉格朗日插值法：
\[A(x)=\sum \limits_{i=1}^n A(x_i) \frac{\prod \limits _{j\ne i}^{n}(x- x_j)}{\prod \limits_{j \ne i}^{n}(x_i-x_j) } \]时间复杂度 \(O(n^2)\)
代码实现：
double res = 0;
for(int i = 1; i <= n; i ++ )
{
double s1 = 1, s2 = 1;
for(int j = 1; j <= n; j ++ )
if(j != i)
{
s1 = s1 * (k - x[j]);
s2 = s2 * (x[i] - x[j]);
}
res = res + y[i] * s1 / s2;
}
复数
基本运算
高中数学基础，不多介绍
\[(a+bi) + (c+di) = (a+c) + (b + d)i \\ (a + bi) - (c + di) = (a - c) + (b - d)i \\ (a + bi)(c + di) = (ac - bd) + (ad + bc)i \\ \frac{a + bi}{c + di} = \frac{ac+bd}{c^2 + d^2} + \frac{bc-ad}{c^2+d^2}i \]欧拉定理
欧拉定理：\(e^{i\theta} = \cos \theta + i\sin \theta\)
根据这个定理，我们可以把复数转化成在复平面上的一堆点来进行一些公式的推导。
算术基本定理
\(n\) 次代数方程 \(f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 = 0(a_n \ne 0 且 n > 0)\) 在复数域上恰好有 \(n\) 个根。
根据这个定理，我们可以找出一个特殊的式子：\(x^n=1\) 的根，在复数域中恰好有 \(n\) 个，设 \(\omega _n = e^{\frac{2\pi}{n}i}\)，则 \(x^n = 1\) 的解集可以记作 \(\{ \omega _n^k|k=0,1,2,\cdots,n-1 \}\)。
单位根
单位根具有一些性质，在 FFT 的时候会有很大用处。
折半引理
\[\omega _{2n}^{k+n} = -\omega _{2n}^{k} \]这个结论在复平面中很好证明，当然也可以用公式来证明：
\[\begin{align} \omega _{2n}^{k+n} &= e^{i\frac{2\pi}{2n}(k+n)}\\ &=e^{i\frac{2\pi}{2n}k+i\pi}\\ &= e^{i\frac{2\pi}{2n}k}\times(-1)\\ &= -\omega _{2n}^k \end{align} \]消去引理
\[\omega _{dn}^{dk} = \omega _{n}^{k} \]证明：
\[\begin{align} \omega _{dn}^{dk} &= e^{i\frac{2\pi}{dn}dk } \\ &=e^{i\frac{2\pi}{n} k} \\ &= \omega ^{k}_{n} \end{align} \]求和引理
\[\sum \limits_{i=0}^{n-1} \omega_n^i= 0 \]证明：
\[\begin{align} \sum \limits_{i=0}^{n-1} \omega_n^i &= \frac{(\omega_n^k)^n -1}{\omega _n^k-1} \\ &=\frac{(\omega _n^n)^k-1}{\omega _n^k-1} \\ &= 0 \end{align} \]FFT
在计算多项式乘法时，系数表示法的复杂度为 \(O(n^2)\)，而点值表示法的复杂度仅为 \(O(n)\)，我们可以尝试将系数表示法转化为点值表示法来进行运算。
我们可以对多项式进行如下变换，将奇数项和偶数项分开：
（前提：\(n\) 为奇数且 \(n+1\) 为 \(2\) 的幂）
\[\begin{align} A(x) &= a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0\\ &= a_nx^n+a_{n-2}x^{n-2}+\cdots+a_1x+a_{n-1}x^{n-1}+a_{n-3}x^{n-3}+\cdots+a_0\\ \end{align} \]将奇数项的系数构成的多项式记为 \(A^{[1]}\)，偶数项记为 \(A^{[0]}\)，就可以得到下面的式子。
\[\begin{align} A(x)&= a_nx^n+a_{n-2}x^{n-2}+\cdots+a_1x+a_{n-1}x^{n-1}+a_{n-3}x^{n-3}+\cdots+a_0\\ &= A^{[0]}(x^2)+xA^{[1]}(x^2) \end{align} \]这时候我们想到，可以对两边的多项式分治处理，但是实数中没有可以把这个平方消掉的数，但是我们可以从复数中找到。
将 \(\omega _{2n}^k\) 带入 \(A\) 中，就可以得到：
\[\begin{align} A(\omega _{2n}^k)&= A^{[0]}(\omega _{2n}^{2k})+\omega _{2n}^kA^{[1]}(\omega _{2n}^{2k}) \\ &= A^{[0]}(\omega _n^k)+\omega _{2n}^kA^{[1]}(\omega _{n}^{k}) \end{align} \]我们就可以对 \(A^{[0]}\) 和 \(A^{[1]}\) 进行递归处理。
同时，在复平面中，\(\omega _{2n}^{n+k}\) 与 \(\omega _{2n}^{k}\) 在多项式 \(A^{[0]}(x^2)\) 和 \(A^{[1]}(x^2)\) 的值相同，因此我们也可以算出该点的值。
\[\begin{align} A(\omega _{2n}^{n+k})&= A^{[0]}(\omega _n^k)+\omega _{2n}^{n+k}A^{[1]}(\omega _{n}^{k})\\ &=A^{[0]}(\omega _n^k)-\omega _{2n}^{k}A^{[1]}(\omega ^k_n) \end{align} \]这就是 FFT 的实现原理。做完 FFT 之后，我们就得到了一个多项式在复数域中 \(2n\) 个点的点值，这时再将两个多项式乘起来，就可以得到乘积的点值表示。
但是如果我们需要的是系数表示，我们则需要做一遍快速傅里叶逆变换，将点值转化为系数。
我们可以列出刚才求点值时的矩阵表示：
\[\begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}=\begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix}\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix} \]我们现在只需要求出逆矩阵即可。由于这个矩阵元素非常特殊（没记错应该叫范德蒙德矩阵），它的逆矩阵即为每一项的倒数再除以变换的长度。
为了取倒数，我们可以做变换：
\[\begin{align} \omega _k^{-1}&= e^{-\frac{2\pi}{k}i } \\ &= \cos (-\frac{2\pi}{k})+\sin(-\frac{2\pi}{k})\\ &=\cos(\frac{2\pi}{k})-\sin(\frac{2\pi}{k}) \end{align} \]我们可以看到与 FFT 时仅差了一个负号，因此我们可以将 FFT 和 IDFT 结合为一题，传参时传入一个 \(opt\)，为 \(1\) 时执行 FFT，为 \(-1\) 时执行 IDFT。
代码实现：
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
const int N = 3e6 + 10;
struct Complex
{
double a, b;
Complex() {a = 0, b = 0;}
Complex(double real, double imag): a(real), b(imag) {}
Complex operator + (const Complex& x) const
{
return Complex(a + x.a, b + x.b);
}
Complex operator - (const Complex& x) const
{
return Complex(a - x.a, b - x.b);
}
Complex operator * (const Complex& x) const
{
return Complex(a * x.a - b * x.b, a * x.b + b * x.a);
}
}F[N], G[N];
void FFT(Complex a[], int lim, int opt)
{
if(lim == 1) return;
Complex a1[lim >> 1], a2[lim >> 1];
for(int i = 0; i < lim; i += 2)
a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
FFT(a1, lim >> 1, opt);
FFT(a2, lim >> 1, opt);
Complex wn = Complex(cos(2.0 * PI / lim), opt * sin(2.0 * PI / lim));
Complex w = Complex(1, 0);
for(int k = 0; k < (lim >> 1); k ++ )
{
a[k] = a1[k] + w * a2[k];
a[k + lim / 2] = a1[k] - w * a2[k];
w = w * wn;
}
}
int n, m;
int main()
{
n = read(), m = read();
for(int i = 0; i <= n; i ++ ) scanf("%lf", &F[i].a);
for(int i = 0; i <= m; i ++ ) scanf("%lf", &G[i].a);
int lim = 1;
while(lim <= n + m) lim <<= 1;
FFT(F, lim, 1), FFT(G, lim, 1);
for(int i = 0; i <= lim; i ++ ) F[i] = F[i] * G[i];
FFT(F, lim, -1);
for(int i = 0; i <= n + m; i ++ ) cout << (int)(F[i].a / lim + 0.5) << ' ';
return 0;
}