Re0 从零开始的多项式之旅
Re0 从零开始的多项式之旅
前の言
温馨提示
- 本文可能不适合 0 基础者观看，因为文章涉及一定的高中数学内容并未特殊说明（主要是懒）；
- 本文可能存在节奏过快的问题，但未提及部分均为读者自推；
- 本文还在更新中；
- 本文可能存在 \(\LaTeX\) 格式书写问题，希望能在评论区中指出批评。
闲话
已经很久没有碰多项式了，一道突然的省选模拟题让死去的回忆复现，这使我充满了决心。当然，可能未曾在此文涉及题目，需要咕咕咕。
快速傅里叶变换
DFT 与 IDFT 的孽缘
多项式算法的最初研究目的就是解决多项式卷积。
但是，直接进行卷积操作不易优化（但可以通过分治——即 Karatsuba 算法，时间复杂度 \(\mathcal O(n^{\log_23})\)），我们想到可以通过一个媒介转化。很显然的，点值表示法是不错的选择，它符合大众的认知（毕竟在初中时期我们就通过待定系数法用三点确定一条二次函数）。
至此，我们就正式引入了 \(\texttt {DFT}\)（离散傅里叶变换，Discrete Fourier Transform）和 \(\texttt {IDFT}\)（离散反傅里叶变换，Inverse Discrete Fourier Transform）。
葬送的单位根
为快速转化，我们将选取一些特殊的点。
warning: 后文将涉及复数和原根内容，如有需要请移步 OI-Wiki
对于 \(\texttt{FFT}\) 而言，它所选取的是 \(w_n^k=e^{\frac{2k\pi}{n}i}\)，其中 \(k\) 从 \(0\) 取到 \(n-1\)，可表示最高次为 \(n-1\) 的多项式。
对于 \(\texttt{NTT}\) 而言，它所选取的是 \(w_n^k=g^{\frac{mod - 1}{n}k}\)，其中 \(g\) 为原根， \(k\) 从 \(0\) 取到 \(n-1\)，可表示最高次为 \(n-1\) 的多项式。
单位根的性质
- \(w_{n}^{x+y}=w_n^xw_n^y\)
- \(w^{2k}_{2n}=w^k_n\)
有了如上基础，就可以来推式子了。
DIT 和DIF 的超距作用
我们令 \(F(k)=\sum\limits_{i=0}^{n-1}a_iw^{ik}_n\)。容易发现，\(F(k)\) 所表示的就是第 \(k+1\) 个我们所选取的点的函数值。
为方便分析，后文的 \(n\) 为 \(2\) 的幂次（实践中也需要先转成 \(2\) 的幂次）。
DIT： 按时域抽取（Decimation-In-Time）
发现，当 \(k < \frac{n}2\) 时，我们令 \(k' = k+\frac n 2\)，则有
显然，由 \(F(k)\) 和 \(F(k')\) 的相似性，我们可以直接递归求解，时间复杂度 \(\mathcal T(n)=2\mathcal T(\frac n 2) + \mathcal O(n)\)，即为 \(\mathcal O(n\log n)\)。
当然，除了递归，我们也可以使用迭代的方式求解。我们需要将 \(a\) 序列变换位置，这被称为 蝴蝶变换（butterfly transform）。
观察发现，在上述式子中 \(a_{2i}\) 和 \(a_{2i+1}\) 分别被移到 \(i\) 和 \(i + \frac{n}{2}\) 位置递归。在二进制上，令 \(2^l=n\)，则 \(a_i\) 最终变换到的位置是 \(i\) 在 \(l\) 位二进制上对称过来的二进制对应的值。例子如下：
13:01101 -> 10110:22
DIF：按频域抽取（Decimation in Frequency）
同理，当 \(k\) 为偶数时，我们令 \(k' = k+1\)，则有
一样的，可以递归，也可以迭代。直接对原序列迭代，最后得到的 \(F\) 是蝴蝶变换后的 \(F\)，再变换一次即可。
\(\texttt{DIT}\) 和 \(\texttt{DIF}\) 的代码之后给出（只有迭代式）。
反方向的钟 —— IDFT 的奇妙冒险
虽然存在一种用矩阵理解的方式，但对于写博客的人来说太折磨了，所以我们用一种简单粗暴的方式来理解并证明它。
事实上，如果我们令 \(\mathrm{DFT}_k(a,n)\) 表示当序列为 \(a\)，长度为 \(n\) 是 \(F(k)\) 的值。同理定义 \(\mathrm{IDFT}_k(a,n)\)，则有
事实上这不好理解，也不易证明，但我还是竭尽码力来书写一番。
证明
即证，\(IDFT_k(\{\sum\limits_{i=0}^{n-1}a_iw^{i}_n,\sum\limits_{i=0}^{n-1}a_iw^{2i}_n,\dots,\sum\limits_{i=0}^{n-1}a_iw^{(n-1)i}_n\},n)=a_k\)。
直接展开，得
显然，当 \(j \not= k\) 时，\(\sum\limits_{i=0}^{n-1}w^{i(j-k)}_n=0\)；而当 \(j=k\) 时，\(\sum\limits_{i=0}^{n-1}w^{i(j-k)}_n=\sum \limits_{i=0}^{n-1}1=n\)。
综上，原式等于 \(a_k\)，即证。
我的代码实现大抵没问题 —— Coding
事实上，当我们从真正意义上将 \(\texttt{DFT}\) 与 \(\texttt{IDFT}\) 联系起来，就可以打出一份完整的代码了：（下为用 \(\texttt{DIT}\) 书写的板子）
// 数组实现版本
using comp = complex<double>;
const double PI = acos(-1);
void FFT(comp *const s, int op) {
for (int i = 0; i < l; ++i) if (i < r[i]) swap(s[i], s[r[i]]); // 做蝴蝶变换
for (int k = 1; k < l; k <<= 1) {
const comp W(cos(PI / k), op * sin(PI / k)); // 单位根
for (int i = 0, K = k << 1; i < l; i += K) {
comp w(1);
for (int j = 0; j < k; ++j, w *= W) {
comp x = s[i + j], y = w * s[i + j + k];
s[i + j] = x + y;
s[i + j + k] = x - y;
}
}
}
if (op == -1) for (int i = 0; i < l; ++i) s[i] /= l;
}
// vector 实现版本（个人推荐）
using comp = complex<double>;
const double PI = acos(-1);
void FFT(vector<comp> &s, int op) {
for (int i = 0; i < l; ++i) if (i < r[i]) swap(s[i], s[r[i]]); // 做蝴蝶变换
for (int k = 1; k < l; k <<= 1) {
const comp W(cos(PI / k), op * sin(PI / k)); // 单位根
for (int i = 0, K = k << 1; i < l; i += K) {
comp w(1);
for (int j = 0; j < k; ++j, w *= W) {
comp x = s[i + j], y = w * s[i + j + k];
s[i + j] = x + y;
s[i + j + k] = x - y;
}
}
}
if (op == -1) for (int i = 0; i < l; ++i) s[i] /= l;
}
然后是蝴蝶变换的递推式（以方便 \(\mathcal O(n)\) 求出变换值）：
for (int i = 0; i < l; ++i)
r[i] = (r[i >> 1] >> 1) | (i & 1 ? l >> 1 : 0);
// 读者易推
浙公网安备 33010602011771号