离散傅里叶变换
实在不知道标题叫啥了，总不能就 DFT 三个字母吧。是最近把 \(\rm DFT\) 的过程自己推了一遍后的浅显理解，以及更快的 \(\rm FFT\) 板子。因为转置原理对于实现的优化实在是太多了，所以我也简单学了一下，不过真的啥都不学不懂。
DFT 在干啥
首先我们要知道 \(\rm DFT\) 是什么，简单来说，它是一组复数 \(\left<a_{0..n-1}\right>\) 到另一组复数 \(\left<b_{0..n-1}\right>\) 的线性变换。式子长这样，其中 \(i\) 是虚数单位：
这个 \(e^{-\frac{i2\pi jk}{n}}\) 其实就是 \(n\) 次单位根的 \(jk\) 次方，\(\omega_n^{jk}\)。可以写出它对应的矩阵：
也叫范德蒙德矩阵。这个矩阵可以在后面的转置原理里用到，这里暂且按下不表。
接下来就来看看它在干什么吧。考虑两个序列 \(\left<f_{0..n-1}\right>,\left<g_{0..n-1}\right>\) 进行 \(\rm DFT\) 后进行点乘，本文约定序列 \(\left<f_{0..n-1}\right>\) 进行 \(\rm DFT\) 后得到的序列为 \(\left<\hat{f}_{0..n-1}\right>\)：
发现也是 \(\rm DFT\) 的形式，考虑设 \(\hat{h}_k=\hat{f}_k\times \hat{g}_k\)，则可以发现：
也就是 \(f_i,g_j\) 在模 \(n\) 意义下的循环卷积。换句话说，我们把 \(\hat{f},\hat{g}\) 点乘就得到了 \(f,g\) 循环卷积后的序列 \(\rm DFT\) 后的结果。如果我们能找到一个逆 \(\rm DFT\)，或者说 \(\rm IDFT\) 的过程，就可以实现循环卷积的计算了。考虑从 \(h_k\) 的表达式入手：
可以发现，\(\rm IDFT\) 的过程既可以通过把 \(\rm DFT\) 的单位根取反得到，也可以通过把序列 \(1\sim n-1\) 处的值反转得到再进行 \(\rm DFT\) 得到。
DFT 该咋做
常见的实现是 \(\rm FFT\)，它的限制是要求循环卷积的模数是 \(2\) 的次幂，通过把这个模数取到大于能得到的最高次幂的方法，我们可以让循环卷积变为普通加法卷积。如果需要实现任意模数的循环卷积，我们需要再推一推式子，即 \(\rm Bluestein's\ algorithm\)。我们分别来看。
首先介绍一下 \(\rm FFT\)。它的基本思路是分治，因为 \(n\) 是 \(2\) 的次幂，所以我们可以每次分成两半来算。考虑设：
从而找到一个分治结构。此外，还可以发现 \(\hat{f}_i\) 与 \(\hat{f}_{i+n/2}\) 式子有一定的相似性：
其实有 \(\hat{f}_{0/1,i+n/2}=\hat{f}_{0/1,i}(i<n/2)\)。从而可以一次计算两个的值，且每次分治长度会减半，从而我们得到一个：
的做法。
递归常数比较大，怎么办呢？模拟一下递归的过程，我们可以找到一个规律。从而可以让我们模拟递归的过程。具体来讲，如果令 \(g_i=f_{\mathrm{rev}_i}\)，其中 \(\mathrm{rev}_i\) 表示将 \(n\) 位二进制数 \(i\) 的位反转得到的数，则第 \(h\) 层就是对：
进行合并。其中 \(0\le i< 2^{h-1}\)。则我们从底向上模拟即可。
实现的时候为了缩小常数我们可以预处理一下单位根。
void FFT(comp* f)
{
for (int i = 0; i < lim; ++i) if (i < rev[i]) std::swap(f[i], f[rev[i]]);
for (int h = 2, t = 1, p = lim / h; h <= lim; h <<= 1, t <<= 1, p >>= 1)
{
for (int j = 0; j < lim; j += h)
for (int q = 0, k = j; q < t; ++q, ++k)
{
comp u = f[k], v = f[k + t] * wn[p * q];
f[k] = u + v; f[k + t] = u - v;
}
}
}
其中 wn[i] 表示 \(\omega_{\rm lim}^i\)。所以 wn[p * q] 为 \(\omega_{\rm lim}^{\rm q\frac{lim}{h}}=\omega_{\rm h}^{\rm q}\)。这里如果是在模 \(p\) 意义下，\(n\) 次单位根就变成了 \(g^{\frac{p-1}{n}}\)，其中 \(g\) 是 \(p\) 的原根。单位根存在，当且仅当 \(n|p-1\)。这也解释了为什么 \(\rm NTT\) 的模数必须满足 \(k2^n+1\) 的形式。
然后我们来看看 \(\rm Bluestein's\ algorithm\)。考虑我们刚刚 \(\rm DFT\) 的过程：
其中利用了 \(ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}\) 这个等式。从而得到了一个差卷积的形式，可以用 \(\rm FFT\) 计算。\(\rm IDFT\) 的过程是类似的。
FFT 实现的优化
利用转置原理对 \(\rm FFT\) 实现进行的常数优化是非常巨大的，且对于减少码量也有帮助。粗体大写表示矩阵，粗体小写表示向量。本部分不会详细系统介绍转置原理。
首先要知道转置定理：
而我们知道，\(\rm DFT\) 的本质是线性变换，所以有：
所以如果我们能找到变换矩阵 \(\bf A\) 的转置，就能以另一种方法计算该线性变换。
对于 \(\rm FFT\)，它的变换矩阵是范德蒙德矩阵，就像我们看到的，它是关于主对角线对称的，所以它的转置等于自身。也就是说转置完等于没转置。看来从整体入手的优化失效了。
考虑把整个大矩阵分成若干小矩阵，即设：
则：
即我们可以把整个过程分成若干小的线性变换，然后倒着执行每个的转置过程。
对应到 \(\rm FFT\) 上，观察我们的过程：
- 对序列做 \(\rm bitrev\) 变换，为模拟递归做准备。
- 模拟递归第 \(n\) 层。
- 模拟递归第 \(n-1\) 层。
- \(\vdots\)
- 模拟递归第 \(1\) 层。
倒过来就是：
- 模拟递归第 \(1\) 层的转置过程。
- \(\vdots\)
- 模拟递归第 \(n\) 层的转置过程。
- 对序列做 \(\rm bitrev\) 变换的转置过程，这个转置过来等于自身。
如果把它和正常的 \(\rm IDFT\) 过程连起来，就能发现 \(\rm bitrev\) 的过程被消掉了。而事实上，由于内存访问不连续等玄学原因，\(\rm bitrev\) 是比较费时间的。从而在这里优化了常数。现在我们需要解决的问题是，模拟递归的转置过程是什么？
考虑模拟递归对应的矩阵。第 \(i\) 层对应的是一个 \(2^{n-i+1}\times 2^{n-i+1}\) 的矩阵。以 \(n-i+1=3\) 为例：
转置！
这样就得到了转置之后的式子。
从而可以写出代码。
inline void NTT(int* f)
{
for (int h = lim, t = lim >> 1, p = 1; h >= 2; h >>= 1, t >>= 1, p <<= 1)
for (int j = 0; j < lim; j += h)
for (int k = j, q = 0; q < t; ++q, ++k)
{
int u = f[k], v = f[k + t];
f[k] = (u + v) % mod; f[k + t] = (ll)w[0][p * q] * (u + mod - v) % mod;
}
}
inline void INTT(int* f)
{
for (int h = 2, t = 1, p = (lim >> 1); h <= lim; h <<= 1, t <<= 1, p >>= 1)
for (int j = 0; j < lim; j += h)
for (int k = j, q = 0; q < t; ++q, ++k)
{
int u = f[k], v = (ll)w[1][p * q] * f[k + t] % mod;
f[k] = (u + v) % mod; f[k + t] = (u + mod - v) % mod;
}
for (int i = 0, inv = ksm(lim, mod - 2); i < lim; ++i) f[i] = (ll)f[i] * inv % mod;
}
其中 w[0/1][i] 分别表示 \(\omega_{\rm lim}^{i}\) 和 \(\omega_{\rm lim}^{-i}\)。真的很快还好写
写完了。我大概是想不开了才会在 \(\rm NOI\) 前总结这东西。
