算法 - FWT & FFT & NTT

反正是用来卷的。

FWT

作用

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

开始

\(C=f(A,B)\)\(A=A_0::A_1\),其中 \(::\) 表示连接, \(A_0,A_1\) 分别表示 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_0::A_1)=fwt(A_0+A_1)::fwt(A_0-A_1)\)

对于 OR,\(fwt(A)=fwt(A_0::A_1)=fwt(A_0)::fwt(A_0+A_1)\)

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

那么 \(C=fwt^{-1}(fwt(A)\circ fwt(B))\),至此 \(fwt\) 完成了。

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

\(n=1\) 时显然成立。

\(n=k\) 时成立,即对于 \(C=f(A,B)\)\(fwt(C)=fwt(A) \circ 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) \circ fwt'(B)\) 的式子,一种方法就是多写一个函数,或者把 A 或 B 取反,那么 g 就满足交换律了。
  • FWT 可以简单地写成非递归的形式。
  • 对于 XOR,有 \(fwt^{-1}(A)=fwt(A)/n\),因为每递归一层就会差 2 倍。
1
2
3
4
5
6
7
8
9
10
11
12
13
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 \in \mathbb R[X]\),求 \(C=A\cdot 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(\omega_n^k)=\sum_n a_i \omega_n^{ki}\),求奇偶项分别求和,设 \(n=2t\),得到 \(F(k)=\sum_t a_{2t}\omega_{t}^{kt}+\omega_n^k \sum_ta_{2t+1}\omega_{t}^{kt}=F_{even}(k)+\omega_n^k F_{odd}(k)\)

对于奇数项和偶数项,分别求出 \(F_{even}\)\(F_{odd}\),但是由于项数只有一半,所以直接获得的只有前 \(t\) 个点。由于 \(t\)\(w_t^k\) 的周期且有 \(\omega_n^k+\omega_n^{k+t}=0\),所以有 \(F(k+t)=F_{even}(k)-\omega_n^k F_{odd}(k)\)

至此 DFT 完成了。

IDFT

但是为了求出 \(C\),还需要一个 DFT 的逆变换 IDFT。

考虑 DFT 的过程,设矩阵 \(M_{ij}=\omega_n^{ij}\),那么 \(DFT(\vec{a})=M \cdot \vec{a}\)

现在需要求 \(IDFT(DFT(\vec{a}))=\vec{a} = M^{-1} \cdot DFT(\vec{a})\)

\(M^{-1}=M'\),那么有 \(M'_{ij}=\frac{1}{n}\omega_n^{-ij}\)

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

至此 FFT 完成了。

实现

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

NTT

与 FFT 唯一的不同是多项式是模意义下的。(\(\mathbb{Z}_m[X]\)

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

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

此处暂略。

后记

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