快速傅里叶变换、数论变换、快速莫比乌斯变换、快速沃尔什变换学习笔记（代码还没补，以及应用）
快速傅里叶变换、数论变换、快速莫比乌斯变换、快速沃尔什变换学习笔记
前置知识
阅读下文需要复数和任意角三角函数等高中数学知识，请先自行学习这些内容。
快速傅里叶变换（FFT）
概述
离散傅里叶变换（Discrete Fourier Transform，DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其 DTFT 的频域采样。快速傅里叶变换（Fast Fourier Transform，FFT）是一种高效实现 DFT 的算法。
DFT 在算法竞赛中常用于加速多项式乘法，并进一步加速多项式其他各种运算来解题。
多项式的系数与点值表示法
系数表示法是用每一项系数来表示多项式。\(f(x)=\sum\limits_{i=0}^{n-1}a_ix^i\iff f(x)=\{a_0,a_1,\cdots,a_{n-1}\}\)。
点值表示法是用 \(n-1\) 次多项式在 \(n\) 个点的取值表示多项式。容易证明 \(n\) 个点的取值可以唯一确定一个不超过 \(n-1\) 次的多项式。\(f(x)=\sum\limits_{i=0}^{n-1}a_ix^i\iff f(x)=\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\}\)。
两个系数表示法的多项式相乘需要计算 \(\mathcal O(n^2)\) 次，而两个点值表示法的多项式相乘只需要在 \(\mathcal O(n)\) 的时间把对应点值相乘即可，于是只要可以在 \(\mathcal o(n^2)\) 的时间实现系数和点值的互转，就可以加速多项式乘法。
离散傅里叶变换（DFT）的作用是将多项式从系数表示法转化为点值表示法，离散傅里叶逆变换（IDFT）的作用是将多项式从点值表示法转化为系数表示法。
了解了 DFT 和 IDFT 的想法，我们发现点值显然不能随便代入，因为这样计算就是 \(\mathcal O(n^2)\) 的了，于是我们希望取一些特殊点的点值。
单位复根
我们想到 \(\pm 1\) 的幂都很好算，就不需要计算 \(x^i\)（\(0\le i < n\)）了。同理，我们发现 \(\pm\textrm{i}\) 的幂也很好算，或者说 \(\omega^k=1\) 中的 \(\omega\) 也可以。
我们定义方程 \(x^n=1\) 在复数意义下的解是 \(n\) 次复根，显然这样的解有 \(n\) 个。我们设 \(\omega_n=\textrm{e}^\frac{2\pi\textrm{i}}{n}=\cos\dfrac{2\pi}{n}+\textrm{i}\sin\dfrac{2\pi}{n}\)，于是 \(x^n=1\) 的解集为 \(\{\omega_n^k\mid k\in[0,n)\cap\Z\}\)。我们称 \(\omega_n\) 为 \(n\) 次单位复根，它是把单位元 \(n\) 等分的第一个角对应的向量。
单位复根有三个重要的性质：
- \(\omega_n^0=\omega_n^n=1\)。
- \(\omega_n^k=\omega_{2n}^{2k}\)。
- \(\omega_{2n}^{k+n}=-\omega_{2n}^k\)。
离散傅里叶变换（DFT）
DFT 的思想是分治地求 \(x=\omega_n^k\) 时 \(f(x)\) 的值，分治的方法是分奇偶次项处理。
注意 DFT 时多项式的次数必须是 \(2^k-1\)，即有 \(2^k\) 项，不足的需要补 \(0\)。
例如对于七次八项式：
进行如下变形：
设函数 \(g,h\)：
则：
代入单位复根：
因此如果我们知道 \(g\left(\omega_{\frac{n}{2}}^k\right)\) 和 \(h\left(\omega_{\frac{n}{2}}^k\right)\)，就可以求出 \(f(\omega_n^k)\) 和 \(f\left(\omega_n^{k+\frac{n}{2}}\right)\)，可以对 \(g,h\) 递归进行 DFT。
时间复杂度为 \(T(n)=2T\left(\dfrac{n}{2}\right)+\mathcal O(n)=\mathcal O(n\log n)\)。
【参考代码先咕着】
位逆序置换
DFT 的分治是递归形式的，时间和空间的常数都较大，FFT 就是通过一些神奇的方式把递归改为递推，使得常数变为原来的一半。
我们对刚刚的七次八项式进行模拟递归：
- \(\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}\)。
- \(\{a_0,a_2,a_4,a_6\},\{a_1,a_3,a_5,a_7\}\)。
- \(\{a_0,a_4\},\{a_2,a_6\},\{a_1,a_5\},\{a_3,a_7\}\)。
- \(\{a_0\},\{a_4\},\{a_2\},\{a_6\},\{a_1\},\{a_5\},\{a_3\},\{a_7\}\)。
那么序列 \(0,4,2,6,1,5,3,7\) 有没有什么特点呢？是有的，把它们写为三位的二进制数，发现就是把二进制高低位翻转了。这一点可以很简单地通过每次考虑最低位的排列方式来证明。这种序列我们称之为位逆序置换，也称蝴蝶变换。
根据位逆序置换的定义，我们显然可以在 \(\mathcal O(n\log n)\) 的时间求出 \(2^k\) 的位逆序置换。
事实上位逆序置换也可以早 \(\mathcal O(n)\) 的时间求出。我们记 \(0\sim 2^k-1\) 的位逆序置换为 \(R(0\cdots 2^k-1)\)，则可以发现：
【参考代码先咕着】
快速傅里叶变换（FFT）
上面说了，FFT 就是 DFT 改成递推版本，我们先对原多项式进行位逆序置换，然后利用倍增合并的方式处理即可。
【参考代码先咕着】
快速傅里叶逆变换（IFFT）
离散傅里叶逆变换（IDFT）
IDFT 可以用 DFT 表示，有两种理解方式。
理解方式一：线性代数
DFT 的过程可以视为将系数序列视为一个列向量，再左乘一个矩阵得到点值列向量：
我们有了点值向量（也就是左侧的列向量），要求系数向量（也就是右侧的列向量），只需要左乘中间矩阵的逆矩阵即可。
中间矩阵十分特殊，它满足：
根据单位复根的性质及欧拉公式：
我们只需要在 DFT 的基础上把单位根 \(\omega_n\) 取成这个数，再把计算结果除以 \(n\) 即可。
理解方式二：单位复根周期性
构造法
原多项式为 \(f(x)=\sum\limits_{i=0}^{n-1}a_ix^i\)，我们已知 \(y_i=f(\omega_n^i)\)（\(0\le i < n\)），要求 \(a_0\cdots a_{n-1}\)。
构造多项式如下：
接下来有两种推导方式，对应两种实现。
推导方式一
设 \(b_i=\omega_n^{-i}\)，则多项式 \(A\) 在 \(x=b_0\cdots b_{n-1}\) 处的点值为 \(A(b_0)\cdots A(b_{n-1})\)。
其中，我们可以把 \(A(b_k)\) 表示为：
我们记 \(S(\omega_n^a)=\sum\limits_{i=0}^{n-1}(\omega_n^a)^i\)。
当 \(a\equiv 0\pmod{n}\) 时，\(S(\omega_n^a)=n\)。
当 \(a\not\equiv 0\pmod{n}\) 时，进行错位相减：
综上，我们知道：
代回原式：
所以我们取单位根为其倒数跑 DFT，再除以 \(n\) 得到的就是原多项式的系数表示。
推导方式二
我们将 \(\omega_n^k\) 代入 \(A\)，与推导方式一类似，最终得到 \(A(\omega_n^k)=\sum\limits_{j=0}^{n-1}a_jS(\omega_n^{j+k})\)，所以 \(A(\omega_n^k)=a_{n-k}\cdot n\)。
所以我们直接对点值做 DFT 后除以 \(n\)，再反转后 \(n-1\) 个元素，得到的也是原多项式的系数表示。
FFT/IFFT 代码实现
代码实现有两种，如下：
【参考代码先咕着】
【参考代码先咕着】
例题
这个我还得多做点，先咕着。
数论变换（NTT）
概述
NTT 是 FFT 在数论基础上的实现，用来解决带模数的多项式乘法。NTT 的运算过程中没有复数和浮点数，常数较小且没有精度误差。
NTT 的详细讲解需要较大量的群论知识，在下文我尝试避开群论知识把 NTT 讲明白。
原根
对于质数 \(p\)，若满足 \(p=q\cdot 2^m+1\)，那么它的原根 \(g\) 满足 \(g^{q\cdot 2^m}\equiv 1\pmod{p}\)。
我们记 \(g_i=g^\frac{p-1}{i}\)，发现 \(g_i\) 与 \(\omega_i\) 有相似的性质。对于正整数 \(n=2^m\) 和整数 \(k\)，有：（下文 \(\implies\) 表示对应关系而非推出）
- \(\omega_n^n=1\implies g_n^n\equiv 1\pmod{p}\)。
- \(\omega_n^k=\omega_{2n}^{2k}\implies g_n^k\equiv g_{2n}^{2k}\pmod{p}\)。
- \(\omega_{2n}^{k+n}=-\omega_{2n}^k\implies g_{2n}^{k+n}\equiv -g_{2n}^k\pmod{p}\)。
- \(\omega_{2n}^n=-1\implies g_{2n}^n\equiv -1\pmod{p}\)。
NTT 模数
常见的 NTT 模数 \(p\) 及其一个原根 \(g\) 有：
- \(p=998244353=7\times 17\times 2^{23}+1,g=3\)。
- \(p=1004535809=479\times 2^{21}+1,g=3\)。
- \(p=469762049=7\times 2^{26}+1,g=3\)。
数论逆变换（INTT）
类似于 IDFT 的推导过程，不过似乎只能用第二种代码实现。
NTT/INTT 代码实现
【参考代码先咕着】
例题
这个我还得多做点，先咕着。
高维前缀和与快速莫比乌斯变换（FMT）
高维前缀和
我们从二维前缀和开始考虑。
二维前缀和有一种很常用的基于容斥的求法：
// 数组下标从 1 开始
for(int i = 1; i <= k; ++i)
for(int j = 1; j <= k; ++j)
a[i][j] = a[i - 1][j] + a[i][j - 1] - a[i - 1][j - 1] + a[i][j];
但这种方法不好推广到更高维，原因是容斥的项数是 \(\mathcal O(2^n)\) 的（其中 \(n\) 为维数），总复杂度就是 \(\mathcal O(k^n2^n)\)。
二维前缀和还有另一种求法，是一维一维求：
// 数组下标从 1 开始
for(int i = 1; i <= k; ++i)
for(int j = 1; j <= k; ++j)
a[i][j] += a[i - 1][j];
for(int i = 1; i <= k; ++i)
for(int j = 1; j <= k; ++j)
a[i][j] += a[i][j - 1];
这种方法是可以推广的。
具体地，对于一个长度为 \(k^n\) 的数组（即：\(n\) 维，每个维度长度为 \(k\)），求它的高维前缀和的方法也是一维一维求：
// 注意：为了方便，与前文不同的是，这里数组下标从 0 开始
// let len := pow(k, n)
for(int j = 0, cur = 1; j < n; ++j, cur *= k)
for(int i = 0; i < len; ++i)
if((i / cur) % k)
a[i] += a[i - cur];
复杂度为 \(\mathcal O(k^nn)\)，可以看到比基于容斥的方法要小很多。正确性的话，可以想象前缀和的求和关系为一张有向无环的网格图，我们在上面按照拓扑序进行 DP，因此是正确的。
一般情况下，每个维度的长度为 \(2\)，高维前缀和也就是子集关系：
// 数组下标从 0 开始
for(int j = 0; j < n; ++j)
for(int i = 0; i < (1 << n); ++i)
if((i >> j) & 1)
a[i] += a[i ^ (1 << j)];
同理，我们还可以写出高维后缀和，也就是超集关系：
// 数组下标从 0 开始
for(int j = 0; j < n; ++j)
for(int i = 0; i < (1 << n); ++i)
if(!((i >> j) & 1))
a[i] += a[i ^ (1 << j)];
Dirichlet 前缀和
狄利克雷前缀和：\(s_k=\sum\limits_{i\mid k}a_i\)。
一种简单的做法是枚举 \(i\)，枚举 \(i\) 的倍数 \(k\)，然后加上去统计答案，这是调和级数的 \(\mathcal O(n\log n)\)。
设 \(k=\prod p_i^{\alpha_i}\)，我们可以把每个质因数 \(p_i\) 视为一维，做高维前缀和，复杂度与埃拉托色尼筛一样为 \(\mathcal O(n\log\log n)\)。
代码：（里面求高维前缀和的函数名叫 fmt，不确定这算不算但不重要）
// Problem: P5495 Dirichlet 前缀和
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P5495
// Memory Limit: 256 MB
// Time Limit: 2000 ms
//
// Powered by CP Editor (https://cpeditor.org)
//By: OIer rui_er
#include <bits/stdc++.h>
#define rep(x,y,z) for(uint x=(y);x<=(z);x++)
#define per(x,y,z) for(uint x=(y);x>=(z);x--)
#define debug printf("Running %s on line %d...\n",__FUNCTION__,__LINE__)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef unsigned int uint;
const uint N = 2e7+5;
uint n, seed, a[N], tab[N], p[N], pcnt, ans;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
uint getnext() {
seed ^= seed << 13;
seed ^= seed >> 17;
seed ^= seed << 5;
return seed;
}
void sieve(uint lim) {
rep(i, 2, lim) {
if(!tab[i]) p[++pcnt] = i;
for(uint j=1;j<=pcnt&&i*p[j]<=lim;j++) {
tab[i*p[j]] = 1;
if(!(i % p[j])) break;
}
}
}
void fmt(uint* a, uint n) {
sieve(n);
rep(i, 1, pcnt) for(uint j=1;j*p[i]<=n;j++) a[j*p[i]] += a[j];
}
int main() {
scanf("%u%u", &n, &seed);
rep(i, 1, n) a[i] = getnext();
fmt(a, n);
rep(i, 1, n) ans ^= a[i];
printf("%u\n", ans);
return 0;
}
概述
FMT 常用于处理子集相关的卷积问题。因为我们可以对集合进行状态压缩，所以 FMT 可以处理按位与、按位或等的卷积。
也就是对于序列 \(a,b\)，求序列 \(c\)，满足：
其中 \(\otimes\) 为 \(\cap\)（交集、按位与）、\(\cup\)（并集、按位或）等。
并集卷积（或卷积）
考虑 FFT 的过程，先将 \(a,b\) 转化为可以直接对应位置相乘的点值表示，再转化回我们想要的系数表示，其中转化的过程用时为 \(\mathcal O(n\log n)\) 低于暴力的 \(\mathcal O(n^2)\)。
那么我们能不能将 \(a,b\) 快速转化为可以直接对应位置相乘的序列 \(\hat{a},\hat{b}\)，求出 \(\hat{c}\) 再快速转化回去呢？
\(a\to\hat{a}\) 的过程被称为快速莫比乌斯变换（FMT），\(\hat{a}\to a\) 的过程被称为快速莫比乌斯逆变换（IFMT）或快速莫比乌斯反演（FMI）。
我们令 \(\hat{a}_i=\sum\limits_{j\subseteq i}a_j\)，则根据 \(c_k=\sum\limits_{i\cup j=k}a_ib_j\)，进行如下推导：
因此上面的 \(\hat{a}\) 可以直接对应位置相乘，符合我们的要求。问题转化为如何 \(a\to\hat{a}\) 以及 \(\hat{a}\to a\)。
先考虑 \(a\to\hat{a}\)。
如果我们暴力枚举 \(\hat{a}\) 定义式 \(\hat{a}_i=\sum\limits_{j\subseteq i}a_j\) 中的 \(i\)，再枚举 \(i\) 的子集 \(j\)，全集大小为 \(\log_2n\)，所以枚举全局的子集的子集的复杂度为 \(\mathcal O(3^{\log_2n})=\mathcal O(n^{\log_23})\approx\mathcal O(n^{1.585})\)，显然不符合我们的要求。
这时候可能会想到，\(\hat{a}_i=\sum\limits_{j\subseteq i}a_j\) 不就是 \(a_j\) 的高维前缀和嘛！于是我们可以 \(\mathcal O(n\log n)\) 求解。
再考虑 \(\hat{a}\to a\)。
我们知道子集反演公式：
于是：
其中 \(|i|\) 表示 \(i\) 集合的大小，状态压缩为二进制后也就是 \(\operatorname{popcnt}(i)\)。
如果你对容斥系数比较敏感的话，会想到把高维前缀和改为“高维前缀差”（瞎起的名字），也就是每次累加贡献改为累加负贡献。
【参考代码先咕着】
交集卷积（与卷积）
类似上面的并集卷积（或卷积），我们构造 \(\hat{a}_i=\sum\limits_{i\subseteq j}a_j\) 即可，也就是把 \(i,j\) 谁是子集谁是超集调换了一下。
代码实现方法和证明方法均同理。
【参考代码先咕着】
\(\gcd\) 卷积和 \(\operatorname{lcm}\) 卷积
也就是 \(c_k=\sum\limits_{\gcd(i,j)=k}a_ib_j\) 和 \(c_k=\sum\limits_{\operatorname{lcm}(i,j)=k}a_ib_j\)。
可以分别构造 \(\hat{a}_i=\sum\limits_{i\mid j}a_j\) 和 \(\hat{a}_i=\sum\limits_{j\mid i}a_j\)，不再赘述。
【参考代码先咕着】
快速沃尔什变换（FWT）
概述
FWT 的思想与 FFT、FMT 相同，也是构造序列 \(a\) 的变换 \(\hat{a}\) 使得可以对应位置相乘。在算法竞赛中，FWT 常用于处理位运算卷积。
也就是对于序列 \(a,b\)，求序列 \(c\)，满足：
其中 \(\otimes\) 为 \(\cap\)（交集、按位与）、\(\cup\)（并集、按位或）等。
并集卷积（或卷积）
我们令 \(\hat{a}_i=\sum\limits_{j\subseteq i}a_j\)，可以用上文 FMT 中的方法证明可以直接对应位置相乘。
先考虑 \(a\to\hat{a}\)。
令 \(a_0,a_1\) 为 \(a\) 的前半段和后半段，也就是 \(a_0\) 是 \(a\) 中下标最高位为 \(0\) 的，\(a_1\) 是下标最高位为 \(1\) 的。
则有：
其中 \(\operatorname{concat}\) 是拼接，\(+\) 是序列对应位置相加，例如 \(\operatorname{concat}([0,1],[2,3])=[0,1,2,3]\)，\([0,1]+[2,3]=[2,4]\)。
于是我们可以分治来解决，而且可以不用递归。
再考虑 \(\hat{a}\to a\)。
由上式可得：
同样分治解决，而且两者还可以合在一起写。
【参考代码先咕着】
交集卷积（与卷积）
类似上面的并集卷积（或卷积），相信大家可以自己推出来了！
【参考代码先咕着】
对称差卷积（异或卷积）
对称差的定义：\(A\Delta B=(A\setminus B)\cup(B\setminus A)=(A\cup B)\setminus(A\cap B)\)。
这应该是算法竞赛中 FWT 最常考的应用（虽说我刚学还没咋见过）。
定义 \(a\odot b=|a\cap b|\bmod 2\)，用状态压缩表示的话就是 \(a\odot b=\operatorname{popcnt}(a\&b)\bmod 2\)。
构造 \(\hat{a}_i=\sum\limits_{i\odot j=0}a_j-\sum\limits_{i\odot j=1}a_j\)，则有：
因此 \(\hat{a}=\operatorname{concat}(\hat{a_0}+\hat{a_1},\hat{a_0}-\hat{a_1})\)，\(a=\operatorname{concat}\left(\dfrac{a_0+a_1}{2},\dfrac{a_0-a_1}{2}\right)\)。
【参考代码先咕着】
再记录一个 \(k\) 进制 FWT 的 博客。
