【学习笔记】FFT
前置知识
重要性质：复数相乘，辐角相加，模长相乘。
即 \(\operatorname{arg}(z_1\times z_2)=\operatorname{arg}(z_1)+\operatorname{arg}(z_2)\)。
那么 \(\operatorname{arg}(z^n)=n\operatorname{arg}(z)\)。
单位根
\(n\) 次幂 为 \(1\) 的复数叫做 \(n\) 次单位根，即方程 \(z^n=1\) 的复数根。\(n\) 次单位根有 \(n\) 个。
因为 \(z^n=1\)，所以 \(\left\vert z\right\vert^n=1\)。而模长一定是个非负实数，所以 \(\left\vert z\right\vert=1\)。也就是说单位根到原点的距离为 \(1\)，一定在单位圆（以原点为圆心，半径为 \(1\) 单位长度的圆）上。
我们把所有 \(n\) 次单位根从 \(1\) 开始，沿逆时针方向记为 \(\omega_n^0,\omega_n^1,\dots\omega_n^{n-1}\)。
由原方程还可以得到 \(\operatorname{arg}((\omega_n^k)^n)=n\operatorname{arg}(\omega_n^k)=\operatorname{arg}(1)\)，所以把 \(1\) 逆时针旋转 \(n\) 次，每次旋转 \(\operatorname{arg}(\omega_n^k)\)，结果还是 \(1\)。那么 \(2\pi\mid n\operatorname{arg}(\omega_n^k)\)。每个单位根又互不相同，所以 \(\operatorname{arg}(\omega_n^k)=\dfrac{k}{n}2\pi\)。它们把单位圆 \(n\) 等分。\(\omega_n^k\) 相当于绕原点把 \(\omega_n^{k-1}\) 逆时针旋转 \(\dfrac{2\pi}{n}\)。
对于 \(k~(k\ge n)\)，定义 \(\omega_n^k=\omega_n^{k\bmod n}\)
性质：
-
\(\omega_n^k=(\omega_n^1)^k\)
-
\(\omega_{nm}^{km}=\omega_n^k\)
因为它们模长都是 \(1\)，辐角也相等。
- \(\omega_n^i\times\omega_n^j=\omega_n^{i+j}\)
用性质 1 容易得证。
- 若 \(n\) 是偶数，\(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)
\(\omega_n^{k+\frac{n}{2}}=\omega_n^k\times\omega_n^\frac{n}{2}\)，\(\operatorname{arg}(\omega_n^\frac{n}{2})=180^\circ\)，相当于把 \(\omega_n^k\) 绕原点旋转 \(180^\circ\)，即关于原点中心对称，实部虚部都取相反数，所以是 \(-\omega_n^k\)。
- \(\sum\limits_{i=0}^{n-1}\omega_n^i=\dfrac{\omega_n^{n-1}\times\omega_n^1-\omega_n^0}{\omega_n^1-1}=0\)
（等比数列求和）
DFT
\(n+1\) 个点确定唯一 \(n\) 次多项式。我们可以用 \(n+1\) 个点表示原多项式，并对点值序列进行一些操作（如乘法），最后由这个点值序列得到一个新的多项式，来代替原多项式之间的操作。
比如 \(F(x)=x^2+2x-1,G(x)=-3x+1\)
在 \(F(x)\) 上取 \((-1,-2)(0,-1)(1,2)(2,7)\)
在 \(G(x)\) 上取 \((-1,4)(0,1)(1,-2)(2,-5)\)
两个点值序列中，横坐标相同的点的纵坐标相乘，得到 \((-1,-8)(0,-1)(1,-4)(2,-35)\)
根据这些点确定一个 \(3\) 次多项式得到 \(-3x^3-5x^2+5x-1\)，恰好等于 \(F(x)\cdot G(x)\)。
DFT 是把系数表达转换成点值表达的过程。我们把单位根代入多项式得到点值序列。
考虑一个 \(n-1\) 次（\(n\) 项）多项式 \(F(x)=\sum\limits_{i=0}^{n-1}F_ix^i\)。
由于 \(n\) 是 \(2\) 的整数幂时比较好做，故 \(n\) 不是 \(2\) 的整数幂时往高次项补 \(0\)。
把它按次数的奇偶性分类：
\(F(x)=F_0+F_2x^2+\dots+F_{n-2}x^{n-2}+F_1x+F_3x^3+\dots+F_{n-1}x^{n-1}\)
令 \(FL(x)=F_0+F_2x+F_4x^2+\dots+F_{n-2}x^{n/2-1},FR(x)=F_1+F_3x+F_5x^2+\dots+F_{n-1}x^{n/2-1}\)
则有 \(F(x)=FL(x^2)+xFR(x^2)\)。
现要把单位根代入 \(F\)。
枚举 \(k~(0\le k<\dfrac{n}{2})\)
\(F(\omega_n^k)=FL((\omega_n^k)^2)+\omega_n^kFR((\omega_n^k)^2)=FL(\omega_{n/2}^k)+\omega_n^kFR(\omega_{n/2}^k)\)
\(F(\omega_n^{n/2+k})=FL((\omega_n^{n/2+k})^2)+\omega_n^{n/2+k}FR((\omega_n^{n/2+k})^2)=FL(\omega_n^{n+2k})+\omega_n^{n/2+k}FR(\omega_n^{n+2k})=FL(\omega_{n/2}^k)-\omega_n^kFR(\omega_{n/2}^k)\)
如果我们知道了 \(x=\omega_{n/2}^k~(0\le k<n/2)\) 时 \(FL(x)\) 和 \(FR(x)\) 的值，就可以算出 \(x=\omega_n^k~(0\le k<n)\) 时 \(F(x)\) 的值。\(FL,FR\) 可以递归求。
求单位根：用三角函数求出 \(\omega_n^1\)，利用性质 1 每次乘上它，防止多次调用三角函数。
\(\large{F(x)=FL(x^2)+xFR(x^2)}\\\large{F(\omega_n^k)=FL(\omega_{n/2}^k)+\omega_n^kFR(\omega_{n/2}^k)}\\\large{F(\omega_n^{n/2+k})=FL(\omega_{n/2}^k)-\omega_n^kFR(\omega_{n/2}^k)}\)
IDFT
有了点值序列还不够，我们还得根据点值求多项式。
令 \(G\) 为点值序列即 \(G_k=F(\omega_n^k)=\sum\limits_{i=0}^{n-1}F_i(\omega_n^k)^i\)
有 \(nF_k=\sum\limits_{i=0}^{n-1}(\omega_n^{-k})^iG_i\)
也就是说，对 \(G\) 再做一遍 DFT，不同的是初始单位根为 \(1\)，每次乘上 \(\omega_n^{-1}\)
而 \(\omega_n^{-1}\) 是 \(\omega_n^{1}\) 的共轭。
不会证
代码实现
根据上文可以写出递归的版本，但常数较大。现在尝试写出非递归版。
考虑 \(0\dots n-1\) 次项系数递归到每一层时的位置。
以 \(n=8\) 为例
0 1 2 3 4 5 6 7
---------------
0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7
写出第一行和最后一行的二进制
| | 原 | 现 |
| 0 | 000 | 000 |
| 1 | 001 | 100 |
| 2 | 010 | 010 |
| 3 | 011 | 110 |
| 4 | 100 | 001 |
| 5 | 101 | 101 |
| 6 | 110 | 011 |
| 7 | 111 | 111 |
我们发现，最后一行 \(i\) 位置上的数其实是位数为 \(\log n\) 时 \(i\) 在二进制下反转得到的数。 于是可以预处理出 \(i\) 在二进制下反转得到的数，也就能得到最后一行。此时他们都是 \(0\) 次，可以直接当做单位根代入得到的值。
然后枚举现在要处理 \(m\) 次多项式（\(m=2,4,8,.\dots,n\)）。把序列分成 \(\dfrac{n}{m}\) 段，枚举段的起始点 \(l\)，在这一段内枚举 \(k~(0\le k<\dfrac{m}{2})\) 计算 \(F(l+k)\) 和 \(F(l+\dfrac{m}{2}+k)\).
代码（多项式乘法）：
const double pi = acos(-1);
int tr[N], n, m;
struct CP {
double x, y;
inline CP operator + (const CP &c) const { return CP(x + c.x, y + c.y); }
inline CP operator - (const CP &c) const { return CP(x - c.x, y - c.y); }
inline CP operator * (const CP &c) const { return CP(x * c.x - y * c.y, x * c.y + y * c.x); }
CP() {}
CP(const double &_x, const double &_y) : x(_x), y(_y) {}
} f[N], g[N];
inline void swap(CP &x, CP &y) {
CP t;
t = x;
x = y;
y = t;
}
void FFT(CP *f, int fl) { // fl=1 DFT ; fl=-1 IDFT
CP w1, buf, tmp;
int now;
rep(i, 0, n - 1) if(i < tr[i]) swap(f[i], f[tr[i]]); // 得到最后一层
for(rei lim = 2; lim <= n; lim <<= 1) { //上文中的 m
now = lim >> 1;
w1 = CP(cos(pi * 2 / lim), fl * sin(pi * 2 / lim)); //w(lim,1)
for(rei i = 0; i < n; i += lim) { //上文中的 l
buf = CP(1, 0);
rep(j, i, i + now - 1) { // 上文中的 l+k
tmp = buf * f[j + now];
f[j + now] = f[j] - tmp;
f[j] = f[j] + tmp;
buf = buf * w1;
}
}
}
}
signed main() {
n = read(); m = read();
rep(i, 0, n) f[i].x = read();
rep(i, 0, m) g[i].x = read();
for(m += n, n = 1; n <= m; n <<= 1) ;
rep(i, 0, n - 1) tr[i] = (tr[i >> 1] >> 1) | (i & 1 ? (n >> 1) : 0); // 预处理二进制反转的结果
FFT(f, 1);
FFT(g, 1);
rep(i, 0, n - 1) f[i] = f[i] * g[i];
FFT(f, -1);
rep(i, 0, m) printf("%d ", (int) (f[i].x / n + 0.49));
putchar(10);
return 0;
}
