浅谈 DFT、IDFT、NTT

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\) 的幂次。

考虑如下一些性质:

  1. \(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\)
  2. \(g_{n}^{0} \equiv g_{n}^{n} \equiv 1\)
    • 显然成立。
  3. \(g_{n}^{l + n} \equiv g_{n}^{l}\)
    • 显然成立。
  4. \(g_{n}^{n/2} \equiv -1\)
    • 证明: 显然 \(g\) 不是 \(\bmod P\) 的二次剩余,则 \(g\)\(\bmod P\) 的非二次剩余。\(\square\)
  5. \((g_{n}^{l})^{k}) \equiv (g_{n}^{k})^{l} \equiv g_{n}^{lk}\)
    • 显然成立。
  6. \(\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;
			}
		}
	}
}
posted @ 2024-10-01 08:17  zsxuan  阅读(6)  评论(0编辑  收藏  举报