算法 - FWT & FFT & NTT

Posted on 2018-08-22 Edited on 2018-08-29 In Algorithms Disqus:

反正是用来卷的。

FWT

作用

用于求 C[k] = ∑_{g(i, j) = k}A[i] ⋅ B[j]，其中 g(i, j) 是一个按位做一个二元运算的函数。

开始

记 C = f(A, B)，A = A₀ : : A₁，其中 : : 表示连接， A₀, A₁ 分别表示 A 的前一半和后一半（首位分别为 0, 1）。

拿 g 为 XOR 为例： $$ C_0 = f(A_0, B_0) + f(A_1, B_1) \\ C_1 = f(A_0, B_1) + f(A_1, B_0) \\ C_0 + C_1 = f(A_0, B_0) + f(A_1, B_1)+f(A_0, B_1) + f(A_1, B_0)=f(A_0+B_0,A_1+B_1) \\ C_0 - C_1 = f(A_0, B_0) + f(A_1, B_1)-f(A_0, B_1) - f(A_1, B_0)=f(A_0-B_0,A_1-B_1) \\ (C_0 + C_1) :: (C_0 - C_1) = f(A_0+B_0,A_1+B_1)::f(A_0-B_0,A_1-B_1) $$ 拿 g 为 OR 为例： $$ C_0 = f(A_0, B_0) \\ C_1 = f(A_0, B_1) + f(A_1, B_0) + f(A_1, B_1) \\ C_0 = f(A_0, B_0) \\ C_0 + C_1 = f(A_0, B_0) + f(A_0, B_1) + f(A_1, B_0) + f(A_1, B_1)=f(A_0+B_0,A_1+B_1) \\ C_0::(C_0 + C_1) = f(A_0, B_0)::f(A_0+B_0,A_1+B_1) $$ 我们递归地定义 fwt 函数，边界为 fwt(a) = a。

对于 XOR，fwt(A) = fwt(A₀ : : A₁) = fwt(A₀ + A₁) : : fwt(A₀ − A₁)。

对于 OR，fwt(A) = fwt(A₀ : : A₁) = fwt(A₀) : : fwt(A₀ + A₁)。

利用之前推出的等式，可以归纳证明 fwt(C) = fwt(A) ∘ fwt(B)，其中 ∘ 为对应位相乘。

那么 C = fwt⁻¹(fwt(A) ∘ fwt(B))，至此 fwt 完成了。

以 XOR 为例简单证明，设 A, B, C 中的元素个数为 n。

n = 1 时显然成立。

若 n = k 时成立，即对于 C = f(A, B) 有 fwt(C) = fwt(A) ∘ fwt(B)。

n = 2k 时 $$ \begin{eqnarray} &&fwt(A)\circ fwt(B) \\ &=&fwt(A_0::A_1)\circ fwt(B_0::B_1) \\ &=&(fwt(A_0+A_1)::fwt(A_0-A_1))\circ(fwt(B_0+B_1)::fwt(B_0-B_1)) \\ &=& (fwt(A_0+A_1) \circ fwt(B_0+B_1)) :: (fwt(A_0-A_1) \circ fwt(B_0-B_1)) \\ &=& fwt(f(A_0+A_1,B_0+B_1))::fwt(f(A_0-A_1,B_0-B_1)) \\ &=& fwt(C_0+C_1)::fwt(C_0-C_1) \\ &=& fwt(C_0::C_1) \\ &=& fwt(C) \end{eqnarray} $$

实现

如果这个运算 g 不满足交换律，你会发现类似 fwt(C) = fwt(A) ∘ fwt^′(B) 的式子，一种方法就是多写一个函数，或者把 A 或 B 取反，那么 g 就满足交换律了。
FWT 可以简单地写成非递归的形式。
对于 XOR，有 fwt⁻¹(A) = fwt(A)/n，因为每递归一层就会差 2 倍。

template<typename T>
void fwt(LL a[], int n, T f) {
    for (int d = 1; d < n; d *= 2)
        for (int i = 0, t = d * 2; i < n; i += t)
            FOR (j, 0, d)
                 f(a[i + j], a[i + j + d]);
}

void XOR (LL& a, LL& b) {
    LL x = a, y = b;
    a = (x + y) % MOD;
    b = (x - y + MOD) % MOD;
}

FFT

作用

已知 A, B ∈ ℝ[X]，求 C = A ⋅ B（多项式乘法）。

设 deg (P) = n，选取 n + 1 个不同的 x_i，得到该多项式的一个点值表示 {(x_i, P(x_i))}，而一个 n + 1 个点值表示能唯一确定一个不超过 n 次的多项式。

如果能够获得 A, B 的点值表示法（对于同样的 x_i），那么可以在线性时间内求出 C 的点值表示法。

而 FFT 正是用于加速这个过程。

开始

DFT

为方便计算，设 $\deg(A)=2^r - 1=N-1, A(x)=\displaystyle \sum_{i=0}^{N-1} a_i x^i$。

选取的 N 个 x_k 为 N 次单位根，即 $\omega_n^k=e^{\frac{2\pi ki}{n}}$。

现在要求 F(k) = A(x_k) = A(ω_n^k) = ∑_na_iω_n^ki，求奇偶项分别求和，设 n = 2t，得到 F(k) = ∑_ta_2tω_t^kt + ω_n^k∑_ta_2t + 1ω_t^kt = F_even(k) + ω_n^kF_odd(k)。

对于奇数项和偶数项，分别求出 F_even 和 F_odd，但是由于项数只有一半，所以直接获得的只有前 t 个点。由于 t 是 w_t^k 的周期且有 ω_n^k + ω_n^k + t = 0，所以有 F(k + t) = F_even(k) − ω_n^kF_odd(k)。

至此 DFT 完成了。

IDFT

但是为了求出 C，还需要一个 DFT 的逆变换 IDFT。

考虑 DFT 的过程，设矩阵 M_ij = ω_n^ij，那么 DFT(a⃗) = M ⋅ a⃗

现在需要求 IDFT(DFT(a⃗)) = a⃗ = M⁻¹ ⋅ DFT(a⃗)。

设 M⁻¹ = M^′，那么有 $M'_{ij}=\frac{1}{n}\omega_n^{-ij}$。

所以只需要对 DFT 算法稍作改动就得到了 IDFT。

至此 FFT 完成了。

实现

DFT 由于要奇偶分别递归，要写成非递归形式可能有些困难。目标是把初始向量重新排列之后变成左右递归而不是奇偶递归，以方便自底向上（其实就是把递归的依据从下标的末尾换成首位），这部分怎么写炫酷就怎么写。

NTT

与 FFT 唯一的不同是多项式是模意义下的。（ℤ_m[X]）

此时需要单位根的替代品，原根，但同时对模数提出了一些要求。

但实现方法和 FFT 是一致的。

此处暂略。

后记

鉴于网上的一些资料不够直观/简洁，于是按自己的思路另外写了一份。