DFT(离散傅里叶变换)
多项式分治。
最早可能是由高斯发现的多项式可以分治,但他的手稿并未作为论文发表。
考虑多项式 \(F(x) = a_0 + a_1 x^{1} + a_2 x^{2} + \cdots + a_{n - 1} x^{n - 1}\) 其中 \(n = 2^{k} \ (k \geq 0)\) 。(任意多项式可以通过高位补 \(0\) 化为这个形式。)
有:
\[\begin{aligned}
F(x) &= a_0 + a_1 x^{1} + a_2 x^{2} + \cdots + a_{n - 1} x^{n - 1} \\
&= a_0 + a_2 x^{2} + a_4 x^{4} + \cdots a_{n - 2} x^{n - 2} + a_1 x + a_3 x^{3} + a_5 x^{5} + \cdots a_{n - 1} x^{n - 1} \\
&= a_0 + a_2 x^{2} + a_4 x^{4} + \cdots a_{n - 2} x^{n - 2} + x (a_1 + a_3 x^{2} + a_5 x^{4} + \cdots a_{n - 1} x^{n - 2}) \\
&= G(x^{2}) + x H(x^{2}) \\
\end{aligned}
\]
最终可以分治为 \(n\) 个关于 \(x^{k}\) 的多项式。
单位根。
(高斯注意到)如果我们将 \(x = e^{i \frac{l}{n} 2 \pi}\) 代入多项式 \(F(x)\) ,则能以 \(O(n \log n)\) 的时间通过分治得到一组关于 \(e^{i \frac{l}{n} 2 \pi}\) 的点值表示。
由欧拉公式 \(e^{i \theta} = \cos \theta + i \sin \theta\) (正确性是你会发现左右两边的泰勒展开一样),将 \(\cos \frac{l}{n} 2 \pi + i \sin \frac{l}{n} 2 \pi\) 的形式写作 \(w_{n}^{l}\) ,即 \(n\) 次单位复根的 \(l\) 次幂。
蝴蝶变换。
(高斯)注意到:
\[\begin{aligned}
F(w_n^{l}) = G(w_n^{2l}) + w_n^{l} H(w_n^{2l}) &= G(w_{n/2}^{l}) + w_n^{l} H(w_{n/2}^{2l}) \\
F(w_n^{l + n/2}) = G(w_n^{2(l + n/2)}) + w_n^{l + n/2} H(w_n^{2(l + n/2)}) = G(w_{n/2}^{l + n/2}) + w_n^{l + n/2} H(w_{n/2}^{l + n/2}) &= G(w_{n/2}^{l}) - w_n^{l} H(w_{n/2}^{l}) \\
\end{aligned}
\]
这展示了 \(F\) 如何由 \(G, H\) 计算。
由于分治最多有 \(\log n\) 层,每层的计算都需要 \(O(n)\) 复杂度。\(DFT\) 的时间是 \(O(n \log n)\) 。
typedef long doube db;
const double PI = acos(-1.);
void DFT(std::vector<CP> *a) {
int n = a.size();
for (int m = 2; m < n; m *= 2) {
db Arg = 1.L / m * 2.L * PI;
std::complex<db> wm(cos(Arg), sin(Arg));
for (int l = 0; l + m - 1 < n; l += m) {
std::complex<db> wmi = wm;
for (int i = l, j = l + m / 2; i < l + m / 2; i++, j++) {
std::complex<db> u = a[i], v = a[j];
a[i] = u + wmi * v;
a[j] = u - wmi * v;
wmi *= wm;
}
}
}
}
IDFT(离散傅里叶逆变换)
注意多项式点值表示的矩阵
\[\begin{aligned}
\left [
\begin{matrix}
1 & x_{0}^{1} & x_{0}^{2} & \cdots & x_{0}^{n - 1} \\
1 & x_{1}^{1} & x_{1}^{2} & \cdots & x_{1}^{n - 1} \\
1 & x_{2}^{2} & x_{2}^{2} & \cdots & x_{2}^{n - 1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n - 1}^{n - 1} & x_{n - 1}^{n - 1} & \cdots & x_{n - 1}^{n - 1} \\
\end{matrix}
\right ]
\cdot
\left [
\begin{matrix}
a_{0} \\
a_{1} \\
a_{2} \\
\vdots\\
a_{n - 1} \\
\end{matrix}
\right ]
=
\left [
\begin{matrix}
F(x_{0}) \\
F(x_{1}) \\
F(x_{2}) \\
\vdots \\
F(x_{n - 1}) \\
\end{matrix}
\right ]
\end{aligned}
\]
\(Vandermonde\) 矩阵可逆,于是:
\[\left [
\begin{matrix}
a_{0} \\
a_{1} \\
a_{2} \\
\vdots\\
a_{n - 1} \\
\end{matrix}
\right ]
=
\begin{aligned}
\left [
\begin{matrix}
1 & x_{0}^{1} & x_{0}^{2} & \cdots & x_{0}^{n - 1} \\
1 & x_{1}^{1} & x_{1}^{2} & \cdots & x_{1}^{n - 1} \\
1 & x_{2}^{2} & x_{2}^{2} & \cdots & x_{2}^{n - 1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n - 1}^{n - 1} & x_{n - 1}^{n - 1} & \cdots & x_{n - 1}^{n - 1} \\
\end{matrix}
\right ]^{-1}
\cdot
\left [
\begin{matrix}
F(x_{0}) \\
F(x_{1}) \\
F(x_{2}) \\
\vdots \\
F(x_{n - 1}) \\
\end{matrix}
\right ]
\end{aligned}
\]
关于范德蒙德矩阵的逆矩阵,有初等爆算和拉格朗日插值法。这里仅说明如何用构造法得到 \(IDFT\) 的逆矩阵。
注意 \(DFT\) 的点值矩阵
\[\begin{aligned}
\left [
\begin{matrix}
1 & (w_{n}^{0})^{1} & (w_{n}^{0})^{2} & \cdots & (w_{n}^{0})^{n - 1} \\
1 & (w_{n}^{1})^{1} & (w_{n}^{1})^{2} & \cdots & (w_{n}^{1})^{n - 1} \\
1 & (w_{n}^{2})^{1} & (w_{n}^{2})^{2} & \cdots & (w_{n}^{2})^{n - 1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & (w_{n}^{n - 1})^{1} & (w_{n}^{n - 1})^{2} & \cdots & (w_{n}^{n - 1})^{n - 1} \\
\end{matrix}
\right ]
\end{aligned}
\]
有
\[\begin{aligned}
&\left [
\begin{matrix}
(w_{n}^{-i})^{0} & (w_{n}^{-i})^{1} & (w_{n}^{-i})^{2} & \cdots & (w_{n}^{-i})^{n - 1}
\end{matrix}
\right ]
\cdot
\left [
\begin{matrix}
(w_{n}^{0})^{j} \\
(w_{n}^{1})^{j} \\
(w_{n}^{2})^{j} \\
\vdots \\
(w_{n}^{n - 1})^{j} \\
\end{matrix}
\right ] \\
&= \sum_{k = 0}^{n - 1} (w_{n}^{-i})^{k} \cdot (w_{n}^{k})^{j} = \sum_{k = 0}^{n - 1} (w_{n}^{-i})^{k} \cdot (w_{n}^{j})^{k} = \sum_{k = 0}^{n - 1} (w_{n}^{j - i})^{k} \\
\end{aligned}
\]
当 \(j = i\) , \(\sum_{k = 0}^{n - 1} 1^{k} = n\) 显然。
当 \(i \neq j\) ,由恒等式 \(x^{n} - 1 = (x - 1)(x^{n - 1} + x^{n - 2} + \cdots + 1)\) ,令 \(x = w_{n}^{l} \ (0 < l < n)\) ,则 \((w_{n}^{l})^{n} = (w_{n}^{n})^{l} = 1\) ,显然 \(w_{n}^{l} \neq 1\) ,则 \(\sum_{k = 0}^{n - 1} (w_{n}^{l})^{k} = 0\) 。
于是
\[\begin{aligned}
\frac{1}{n}
\cdot
\left [
\begin{matrix}
1 & (w_{n}^{-0})^{1} & (w_{n}^{-0})^{2} & \cdots & (w_{n}^{-0})^{n - 1} \\
1 & (w_{n}^{-1})^{1} & (w_{n}^{-1})^{2} & \cdots & (w_{n}^{-1})^{n - 1} \\
1 & (w_{n}^{-2})^{1} & (w_{n}^{-2})^{2} & \cdots & (w_{n}^{-2})^{n - 1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & (w_{n}^{-(n - 1)})^{1} & (w_{n}^{-(n - 1)})^{2} & \cdots & (w_{n}^{-(n - 1)})^{n - 1} \\
\end{matrix}
\right ]
\cdot
\left [
\begin{matrix}
1 & (w_{n}^{0})^{1} & (w_{n}^{0})^{2} & \cdots & (w_{n}^{0})^{n - 1} \\
1 & (w_{n}^{1})^{1} & (w_{n}^{1})^{2} & \cdots & (w_{n}^{1})^{n - 1} \\
1 & (w_{n}^{2})^{1} & (w_{n}^{2})^{2} & \cdots & (w_{n}^{2})^{n - 1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & (w_{n}^{n - 1})^{1} & (w_{n}^{n - 1})^{2} & \cdots & (w_{n}^{n - 1})^{n - 1} \\
\end{matrix}
\right ]
=
\left [
\begin{matrix}
1 & 0 & 0 & \cdots & 0 \\
0 & 1 & 0 & \cdots & 0 \\
0 & 0 & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & 1 \\
\end{matrix}
\right ]
\end{aligned}
\]
于是得到了用于 \(IDFT\) 的逆矩阵。
注意 \(w_{n}^{-l} = \overline{w_{n}^{l}}\) 。于是 \(IDFT\) 时,只需将用共轭单位根进行 \(DFT\) 能得到答案的 \(n\) 倍。只需最后让答案除以 \(n\) 。
\(DFT + IDFT\) 即数论函数意义下的 \(FFT\) 。
FFT 精度问题。
显然浮点单位根幂乘上去精度很可能爆炸,最好是使用 \(long\ double\) 的浮点数。
NTT
数论函数中,当所需要的结果很大。无法使用普通数据类型存储,需要对结果取模。
显然不能用 \(FFT\) 直接做乘法,至少单位根累乘的过程中无法在浮点数上取模。
特殊模数模意义下的 \(FFT\) 又叫 \(NTT\) ,全程快速数论变换。
考虑一个质数 \(P\) 同时也为 \(2\) 的幂次。令 \(\mod p\) 的一个原根为 \(g\) ,钦定 \(g^{\frac{P - 1}{n} \times l} = g_{n}^{l}\) 。考虑 \(n\) 为 \(2\) 的幂次。
考虑如下一些性质:
- \(g_{n}^{2l} \equiv g_{n / 2}^{l}\) 。
- 证明: \(g_{n}^{2l} = g^{\frac{P - 1}{n} 2l} = g^{\frac{P - 1}{n / 2} l}\) 。由 \(n\) 为 \(2\) 的幂次显然成立。\(\square\)
- \(g_{n}^{0} \equiv g_{n}^{n} \equiv 1\) 。
- \(g_{n}^{l + n} \equiv g_{n}^{l}\) 。
- \(g_{n}^{n/2} \equiv -1\) 。
- 证明: 显然 \(g\) 不是 \(\bmod P\) 的二次剩余,则 \(g\) 是 \(\bmod P\) 的非二次剩余。\(\square\)
- \((g_{n}^{l})^{k}) \equiv (g_{n}^{k})^{l} \equiv g_{n}^{lk}\) 。
- \(\sum_{k = 0}^{n - 1} g_{n}^{l} \equiv 0 \ s.t. 0 < l < n\) 。
- 证明: 考虑 \(x^{n} - 1 = (x - 1)(x^{n - 1} + x^{n - 2} + \cdots + 1)\) 。\((g_{n}^{l})^{n} = 1\) 显然,\(g_{n}^{l} \neq 1\) 显然,于是 \(\sum_{k = 0}^{n - 1} (g_{n}^{l})^{k} = 0\) 。\(\square\)
于是 \(g^{\frac{P - 1}{n} l} = g_{n}^{l}\) 具备了 \(w_{n}^{l}\) 在 \(FFT\) 中的所有性质。
void NDFT(std::vector<i64> *a) {
int n = a.size();
for (int m = 2; m < n; m *= 2) {
i64 gm = g; // g^{-1} = ksm(g, P - 2);
for (int l = 0; l + m - 1 < n; l += m) {
i64 gmi = g;
for (int i = l, j = l + m / 2; i < l + m / 2; i++, j++) {
i64 u = a[i], v = a[j];
a[i] = (u + gmi * v) % P;
a[j] = ((u - gmi * v) % P + P) % P;
gmi *= g;
}
}
}
}