多项式乘法（FFT）
多项式乘法（FFT）
前言
先说几个比较容易混淆的概念的缩写吧
DFT：离散傅里叶变换 \(O(n^{2})\) 计算多项式乘法
FFT：快速傅里叶变换 \(O(n\log(n))\) 计算多项式乘法
FNTT/NTT：快速傅里叶变换的优化版 优化常数及误差
FWT：快速沃尔什变换 利用类似 FFT 的东西解决一类卷积问题
MTT：毛爷爷的 FFT 非常 nb/任意模数
FMT：快速莫比乌斯变化
多项式
称一个关于 \(x\) 的式子
为一个 \(n\) 次多项式，其中 \(a_{i}\) 为常数。称 \(n\) 为 \(f(x)\) 的次数。
对于一个 \(n\) 次函数的解析式，则需要 \(n+1\) 个点的坐标。这是因为一个 \(n\) 此函数一共有 \(n+1\) 个系数。
也就是说，只要给定了 \(n+1\) 个点的坐标，就可以唯一确定一个 \(n\) 此函数的解析式，进而也就可以确定这个多项式。
那么给出 \(n+1\) 个互相不相同的点的坐标来确定一个多项式的形式，称为多项式的 点值表示法，对应的，写为 \(f(x) =\sum\limits_{i = 0}^{n}a_{i}\times x^{i}\) 形式的多项式被称为多项式的系数表示法
复数
关于复数基本定义不再赘述。
辐角：假设以逆时针为正方向，从 \(x\) 轴正半轴到已知向量的转角的有向角叫做辐角
在复数的几何意义中，复数相乘就是辐角相加，模长相乘
在此，我们会用到欧拉公式：\(re^{i\theta} = a+bi\)，其中 \(\theta = \arctan(\frac{b}{a}), r = \sqrt{a^2+b^2}\)
在这里，我们用它的一种具体形式：\(e^{i\theta} = \cos x + i\sin x\)
单位根
在复平面上，以原点为圆心，1 为半径做圆，所得的圆叫做单位圆。以原点为起点，圆的 \(n\) 等分点为终点，做 \(n\) 个向量，设辐角为正且最小的向量对应的复数为 \(\omega _{n}\) 称为 \(n\) 次单位根。
根据复数乘法的运算法则，其余 \(n - 1\) 个复数为 $\omega_{n}{2}, \omega_{n}, \ \omega_{n}{4},...,\omega_{n} $
所以由欧拉公式，\(\omega_{n}^{k} = \cos(\frac{2\pi}{n} \times k) + i \sin (\frac{2\pi}{n} \times k) = e^{i \frac{2\pi}{n} \times k}\)
在代数中，如果 \(z^{n} = 1\) ，我们把 \(z\) 称为 \(n\) 次单位根
单位根的性质
-
\(\omega_{2n}^{2k} = \omega_{n}^{k}\)
由单位根定义就可以得到，也可以根据欧拉公式来证明 -
设 \(n\) 是一个偶数，且 \(n = 2m\) ，有：
\[(\omega_{n}^{k})^{2} = \omega_{m}^{k} \\ \omega_{n}^{m + k} = -\omega_{n}^{k} \\ \]
一式的证明，可以考虑两个 \(\omega_{n}^{k}\) 相乘，辐角相加，为之前辐角的两倍，模长为 \(1\times 1 = 1\)
二式的证明，左边可以写成 \(\omega_{n}^{k} \times \omega_{n}^{m}\)，即将 \(\omega_{n}^{k}\) 绕原点旋转 \(\pi\) 弧度，所以横纵坐标都是原来的相反数。
算法
如果对两个多项式做乘法，运用系数表示法，需要 \(O(n^{2})\) 的时间复杂度，而如果已知两个多项式的点值表示法，则只需要 \(O(n)\) 的时间复杂度，因为只需要将对应点的纵坐标相乘就可以了。
但是我们将一个多项式从系数表示法转变为点值表示法（称为求值）需要 \(O(n^{2})\) 的时间复杂度（因为每一个横坐标都需要 \(O(n)\) 的时间去计算纵坐标），而将一个点值表示法改为系数表示法（称为插值）则需要 \(O(n^{3})\) 的复杂度来做高斯消元。
所以我们加速多项式乘法的过程，变为了：
系数表示 $\Rightarrow $ 点值表示 \(\Rightarrow\) 乘法 \(\Rightarrow\) 点值表示 \(\Rightarrow\) 系数表示
而傅里叶变换就是从系数表示转换到点值表示的桥梁，我们观察到，点值表示下所取的点是可以任取的，所以我们想利用这一个性质，利用一些特殊值，看是否可以加速运算。
求出一个 \(n\) 次多项式在每个 \(n\) 次单位根下的点值的过程，被称为 离散傅里叶变换（Discrete Fourier Transform, DFT），而将这些点值重新插值成系数的过程，叫做离散傅里叶逆变换(Inverse Discrete Fourier Transform, IDFT)
我们下面进行的变换的多项式为 \(n - 1\) 次多项式 \(A(x) =\sum\limits_{i = 0}^{n - 1}a_i\times x^{i}\) 。其中，\(n\) 是 2 的整数次幂。如果不是的话，就向 \(A(x)\) 的更高次数位 \(n\) 补充 \(a_n = 0\)，令其变为 \(n\) 次多项式，一直进行到其次数 \(+1\) 的值是 2 的整数幂，取 \(n\) 等于其次数，\(m = \frac{n}{2}\)
DFT
考虑求出一个长度为 \(n\) 的数列 \({b_{i}}\)，这个数列的第 \(k\) 项的为 \(A(x)\) 在 \(n\) 次单位根的 \(k\) 次幂处的点值，因此有：
这个过程是 \(O(n^{2})\) 的，我们考虑用快速傅里叶变换(Fast Fourier Transform, FFT) 来优化这个过程。
FFT
我们考虑对于 \(A(x)\) 按照系数角标的奇偶性分类，即:
对于上式的后半部分，我们考虑提出一个 \(x\)，得到：
设前面的多项式是 \(A_{0}(x)\)，后面的是 \(A_1(x)\)，那么可以得到
如果求出了 \(A_0\) 和 \(A_1\) 在各点的点值，那么可以根据上式 \(O(1)\) 计算 \(A(x)\) 在各个点的点值了。而求 \(A_0\) 和 \(A_1\) 的过程和求 \(A\) 的过程完全一致，因此可以递归处理。但是很遗憾，\(A_0\) 和 \(A_1\) 各需要递归一次，根据主定理，复杂度是 \(O(n^{2})\) 的
但是考虑我们求的是在 \(n\) 次单位根的各个幂下的点值，在单位根性质里面的第二条中，使得我们考虑小于 \(m\) 次的点值和大于 \(m\) 次的点至之间的关系
对于 \(0 \leq k < m\) ，我们有：
根据上面的第一个公式化简得到：
这只是前半部分的次数，我们考虑后半部分的次数：
根据第二个公式和第一个公式，得到：
我们惊喜地发现，大于 \(m\) 次的点值可以由 \(A_0\) 和 \(A_1\) 在小于 \(m\) 次的点值求出。只要求出了 \(A_0\) 和 \(A_1\) 在小于 \(m\) 次的点值，就可以线性求出 \(A\) 在整个 \(n\) 次幂处的点值。而求 \(A_0\) 和 \(A_1\) 在小于 \(m\) 次的点值也是可以递归求解的。
考虑时间复杂度：递推关系为 \(T(n) = 2T(n/2) + O(n)\) ，所以 \(T(n) = n\log (n)\)
由此我们得到了快速计算 DFT 的方法，这种方法就叫做 FFT
上面推导出 \(A(\omega_{n}^{k})\) 和 \(A(\omega_{n}^{m+k})\) 的值的式子被称为蝴蝶操作(Butterflt Operation) 。这个名字的由来是如果将 \(A\) 在各点的点值画成一个长条形的数组，在数组下面依次是 \(n\) 本愿单位根的 \(0到 m - 1\) 次幂，求值过程可以画成 \(\omega_{n}^{0}\) 连一条边向 \(A(\omega_{n}^{0})\) ，\(A(\omega_{n}^{\frac{n}{2} - 1})\) ，而 \(\omega_{n}^{1}\) 连向 \(A(\omega_{n}^{1}), A(\omega_{n}^{\frac{n}{2}})\) ，以此类推，画出的图形形状如同蝴蝶的翅膀。
IDFT
至此，我们已经有了 \(O(n \log n)\) 的算法计算两个系数表示法的多项式相乘后的点值表示。接下来我们只需要用 \(O(n \log n)\) 的时间复杂度完成插值的过程，就可以得到一个完整的 \(O(n \log n)\) 的系数型多项式乘法计算了。
我们目前已知一个 \((n - 1)\) 次多项式 \(A(x) = \sum\limits_{i = 0}^{n - 1} a_ix^{i}\) 进行了离散傅里叶变换后的点值 \({b_i}\) ，即
现在试图还原系数数列 \({a_i}\)
推导过程比较复杂，我们直接给出结论：
下面证明上面这个式子是 DFT 式子(即上面计算的 \(b_{k}\) 的式子) 逆变换
我们首先将 DFT 的式子带入 \(\sum\limits_{i = 0}^{n - 1} b_i\omega_{n}^{-ki}\)
得到上式为：
考虑分类讨论
-
当 \(j = k\) 时：
则 \(j - k = 0\) ，因此 \(\omega_{n}^{i(j - k)} = 1\) ，于是\[\sum\limits_{i = 0}^{n - 1}\omega_{n}^{i(j - k)} = \sum\limits_{i = 0}^{n - 1} 1 = n \] -
当 \(j \neq k\) 时：
显然 \(\abs{j - k} < n\)，因此 \(\omega_{n}^{j - k} \neq 1\)，所以 \(\sum\limits_{i = 0}^{n - 1}\omega_{n}^{(j - k)i}\) 是一个公比不为 1 的定比数列前缀和，由等比数列求和公式可以得到：\[\sum\limits_{i = 0}^{n - 1}\omega_{n}^{(j - k)i} = \frac{1 - (\omega_{n}^{j - k})^{n}}{1 - \omega_{n}^{j - k}} \]根据复数的几何性质，易证对于所有的 \(n\) 次单位根 \(T\)，\(T^{x+n} = T^{x}\)，其中 \(x\) 为任意整数
\(\sum\limits_{i = 0}^{n - 1} \omega_{n}^{(j - k)i} = \frac{1 - (\omega_{n}^{j - k})^{n}}{1 - \omega_{n}^{j - k}} = \frac{1 - 1}{1 - \omega_{n}^{j - k}} = 0\)
第二中情况的证明被称为 消去引理
综上讨论，当 \(j \neq k\) 时，因为另一个因数是 0，前面的 \(\sum a_{j}\) 不会对式子产生贡献，当$ j = k$ 的时候，会对答案产生 \(n\) 倍的贡献。所以原式可以写为：
类似地，可以证明将 IDFT 的式子带入 DFT 也是可以使等式成立的。
所以我们证明了变换：
是 DFT 的逆变换，称为 IDFT。我们可以通过这个式子求出这个多项式的系数表示法。
IFFT
下面的问题是如何用较低的复杂度计算 \(B(x) = \sum\limits_{i = 0}^{n - 1} b_{i} \times x^{i}\) ，在 \(\omega_{n}^{-ki}\) ，其中 \(0 \leq k < n\) 处的点值。
而 \(\omega_{n}^{-k}\) 可以看做 \(n\) 次本原单位根每次逆时针旋转本原单位根角的弧度，因此 \(\omega_{n}^{-k}\) 和 \(\omega_{n}^{k}\) 是一一对应的。具体的，\(\omega_{n}^{-k} = \omega_{n}^{k + n}\) 。因此，我们只需要使用 FFT 的方法，求出 \(B(x)\) 在 \(\omega_{n}\) 各个幂次下的值，然后数组反过来，即令 \(a_{k} = \frac{1}{n} \sum\limits_{i = 0}^{n} B(\omega_{n}^{n - k})\) 即可。
这一步快速计算插值的过程叫做 快速傅里叶逆变换(Inverse Fast Fourier Transform , IFFT)
至此，我们得到了一个时间复杂度为 \(O(n \log n)\) 的多项式乘法计算方法。
Code
void FFT(std::complex<double> *A, int N) {
if (N == 1) {
return;
}
int M = N >> 1;
std::complex<double> A0[M], A1[M];
for (int i = 0; i < M; ++i) {
A0[i] = A[i << 1];
A1[i] = A[(i << 1) | 1];
}
FFT(A0, M); FFT(A1, M);
auto W = std::complex<double>(cos(1.0 * PI / M), sin(1.0 * PI / M)), w = std::complex<double>(1.0, 0.0);
for (int i = 0; i < M; ++i) {
A[i] = A0[i] + w * A1[i];
A[i + M] = A0[i] - w * A1[i];
w *= W;
}
}
以及 IFFT：
void IFFT(std::complex<double> *A, int N) {
FFT(A, n)
std::reverse(A + 1, A + N);
}
需要注意的是，因为 \(\omega_{n}^{0} = \omega_{n}^{n}\) ，所以是第 1 个点值和第 \(N-1\) 个交换，而不是第 0 个点值和第 \(N - 1\) 个交换
参考
picks.logdown.com/posts/177631-fast-fourier-transform
