快速傅里叶变换
模板题：【模板】多项式乘法（FFT）
对于两个多项式 \(f,g\)，\(f\) 的项数为 \(n\)，\(g\) 的项数为 \(m\)，FFT 可以 \(O((n+m)\log(n+m))\) 地求它们的卷积。
多项式点值表示法
引理：给定 \(n\) 个点值 \(f(x_0),f(x_1),f(x_2),\dots,f(x_{n-1})\)（\(x_i\) 互不相同），可确定唯一一个多项式 \(f(x)=a_0+a_1x+a_2x^2+a_3x^3+\dots+a_{n-1}x^{n-1}\)。
证明：
对于给定 \(n\) 个点值 \(f(x_0),f(x_1),f(x_2),\dots,f(x_{n-1})\)，可得方程：
\[\begin{cases} a_0+a_1x_0+a_2x_0^2+a_3x_0^3+\dots+a_{n-1}x_0^{n-1}=f(x_0)\\ a_0+a_1x_1+a_2x_1^2+a_3x_1^3+\dots+a_{n-1}x_1^{n-1}=f(x_1)\\ a_0+a_1x_2+a_2x_2^2+a_3x_2^3+\dots+a_{n-1}x_2^{n-1}=f(x_2)\\ \dots\\ a_0+a_1x_{n-1}+a_2x_{n-1}^2+a_3x_{n-1}^3+\dots+a_{n-1}x_{n-1}^{n-1}=f(x_{n-1}) \end{cases}\]这是一个由 \(n\) 个本质不同的 \(n-1\) 元 \(n-1\) 次方程组成的方程组，有唯一解。
多项式转点值
考虑用 \(n\) 个点值 \(f(x_0),f(x_1),f(x_2),\dots,f(x_{n-1})\) 表示多项式 \(f\)，现在需要快速求出这 \(n\) 个点值。
假设 \(f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7\)，有：
\[\begin{aligned} f(x)&=(a_0+a_2x^2+a_4x^4+a_6x^6)+(a_1x+a_3x^3+a_5x^5+a_7x^7)\\ &=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6) \end{aligned}\]令 \(f_1(x)=a_0+a_2x+a_4x^2+a_6x^3,f_2(x)=a_1+a_3x+a_5x^2+a_7x^3\)，有 \(f(x)=f_1(x^2)+x\times f_2(x^2)\)。
上式仍只能每次求解一个 \(f(x_i)\)，考虑用单位根作为点值，有：
\[\begin{aligned} &\begin{aligned} f(\omega_n^k)&=f_1((\omega_n^k)^2)+\omega_n^k\times f_2((\omega_n^k)^2)\\ &=f_1(\omega_n^{2k})+\omega_n^k\times f_2(\omega_n^{2k})\\ &=f_1(\omega_{\frac{n}{2}}^k)+\omega_n^k\times f_2(\omega_{\frac{n}{2}}^k) \end{aligned}\\ &\begin{aligned} f(\omega_n^{k+\frac{n}{2}})&=f_1((\omega_n^{k+\frac{n}{2}})^2)+\omega_n^{k+\frac{n}{2}}\times f_2((\omega_n^{k+\frac{n}{2}})^2)\\ &=f_1((-\omega_n^k)^2)-\omega_n^k\times f_2((-\omega_n^k)^2)\\ &=f_1(\omega_n^{2k})-\omega_n^k\times f_2(\omega_n^{2k})\\ &=f_1(\omega_{\frac{n}{2}}^k)-\omega_n^k\times f_2(\omega_{\frac{n}{2}}^k) \end{aligned} \end{aligned} \]这样就把所求的点值联系了起来，可以递归求解：
\(f(\omega_1^0)\to f(\omega_2^{0\sim 1})\to\dots\to f(\omega_{\frac{n}{4}}^{{0\sim\frac{n}{4}-1}})\to f(\omega_{\frac{n}{2}}^{{0\sim\frac{n}{2}-1}})\to f(\omega_n^{{0\sim n-1}})\)
时间复杂度 \(O(n\log n)\)。
该算法要求 \(n\) 为 \(2\) 的幂次，对于不是 \(2\) 的幂次的情况需要找到一个 \(2^k\ge n\)，把 \(n+1\sim 2^k\) 的系数补为 \(0\)。
但递归常数巨大，不实用。
考虑先求出底层的 \(f(\omega_1^0)\)，再一层一层向上计算，就省去了递归的大常数。
于是需要求出底层 \(a\) 的排列。
假设 \(f\) 的各项系数为 \(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\)，递归法会这样划分：
-
\((a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)\)
-
\((a_0,a_2,a_4,a_6),(a_1,a_3,a_5,a_7)\)
-
\((a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)\)
-
\((a_0),(a_4),(a_2),(a_6),(a_1),(a_5),(a_3),(a_7)\)
假设一开始每个数都在位置 \(0\)，若 \(i\) 的二进制下第 \(i\) 位为 \(1\)（所有位为 \(1\sim k\)），那么 \(i\) 在排列中会被往后移 \(2^{k-i}\) 位。
因此，在二进制下将 \(i\) 翻转，就得到了它在底层排列的位置。
点值转多项式
令 \(F(x)=f(x_0)+f(x_1)x+f(x_2)x^2+\dots+f(x_{n-1})x^{n-1}\)，有：
\[\begin{aligned} F(\omega_n^{-k})&=\sum_{i=0}^{n-1}f(\omega_n^i)(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\Big(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j\Big)(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\Big(\sum_{j=0}^{n-1}a_j\omega_n^{ij}\Big)\omega_n^{-ik}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j-k)}\\ &=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i \end{aligned}\]令 \(g(\omega_n^t)=\displaystyle\sum_{i=0}^{n-1}(\omega_n^t)^i\)
分类讨论 \(g(\omega_n^t)\) 的值：
-
\(t=0\)
\(g(\omega_n^0)=\displaystyle\sum_{i=0}^{n-1}(\omega_n^0)^i=\sum_{i=0}^{n-1}1^i=n\)
-
\(t\ne 0\)
根据等比数列求和公式，\(g(\omega_n^t)=\dfrac{(\omega_n^t)^n-1}{\omega_n^t-1}=\dfrac{1-1}{\omega_n^t-1}=0\)
故有：
\[g(\omega_n^t)=\begin{cases} n,&\text{}t=0\\ 0,&\text{}t\ne 0 \end{cases}\]带入原式：
\[\begin{aligned} F(\omega_n^{-k})&=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ &=\sum_{j=0}^{n-1}a_jg(\omega_n^{j-k})\\ &=a_k\times n \end{aligned}\]也就是说 \(F(\omega_n^{-k})=\dfrac{f(\omega_n^k)}{n}\)。
因此，将点值作为 \(F\) 的系数，再求一遍 \(F\) 的点值，除以 \(n\)，就得到了 \(f\) 的各项系数。
多项式运算转点值运算
显然，两个多项式运算可以转化为点值运算。
例如两个多项式相乘（卷积），就相当于它们的对应点值相乘。其他运算同理。
因此，利用 FFT 将难以直接运算的多项式转为极好运算的点值，就可以方便地进行多项式运算。
代码
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
int n, m, rev[4000005];
struct cpx {
double x, c;
};
cpx operator+(cpx x1, cpx x2) {
return {x1.x + x2.x, x1.c + x2.c};
}
cpx operator-(cpx x1, cpx x2) {
return {x1.x - x2.x, x1.c - x2.c};
}
cpx operator*(cpx x1, cpx x2) {
return {x1.x *x2.x - x1.c * x2.c, x1.x *x2.c + x2.x * x1.c};
}
void change(int len, cpx *f) {
for (int i = 0; i < len; i++) {
rev[i] = rev[i >> 1] >> 1;
if (i & 1)
rev[i] |= len >> 1;
}
for (int i = 0; i < len; i++)
if (i < rev[i])
swap(f[i], f[rev[i]]);
}
void FFT(int op, int len, cpx *f) {
cpx w, tw, t1, t2;
change(len, f);
for (int t = 2; t <= len; t <<= 1) {
if (op)
tw = {cos(PI * 2 / t), sin(PI * 2 / t)};
else
tw = {cos(PI * 2 / t), sin(-PI * 2 / t)};
for (int i = 0; i < len; i += t) {
w = {1, 0};
for (int j = i; j < i + t / 2; j++) {
t1 = f[j], t2 = f[j + t / 2];
f[j] = t1 + w * t2, f[j + t / 2] = t1 - w * t2;
w = w * tw;
}
}
}
if (!op) {
for (int i = 0; i < len; i++)
f[i].x /= len;
}
}
cpx f[4000005], g[4000005];
int main() {
cin >> n >> m;
for (int i = 0; i <= n; i++)
cin >> f[i].x;
for (int i = 0; i <= m; i++)
cin >> g[i].x;
int len = 1;
while (len < n + m + 1)
len <<= 1;
FFT(1, len, f);
FFT(1, len, g);
for (int i = 0; i < len; i++)
f[i] = f[i] * g[i];
FFT(0, len, f);
for (int i = 0; i < n + m + 1; i++)
cout << (int)(f[i].x + 0.5) << ' ';
//四舍五入，避免精度误差
return 0;
}