『学习笔记』FFT与NTT
作者：互联网
这个傻逼目录为什么挂了啊？？
中文名：快速傅里叶离散变换、快速数论变换
英文名：fast fast tle FFT,NTT
之前看这东西看了好久没看懂，回来看原来只是几个无脑操作
前置芝士：三角函数基础(FFT)、原根(NTT)
\(\mathtt {FFT}\) 部分
1. 胡扯
假如现在给你两个多项式 \(f,g\) ，怎么算出它们的加法卷积 \(c\) 呢？
由于 \(c[k]=\sum\limits_{i+j=k}^{}f[i]g[j]\) ，暴力算这玩意显然是 \(O(n^2)\) 的
于是我们需要快速TLE快速傅里叶离散变换来帮我们解决这个问题
2. 继续胡扯
\(\mathtt{DFT\ \&\ IDFT}\)
小学二年级的时候，数学老师告诉我们，在平面直角坐标系上 \(n+1\) 个点可以确定一个 \(n\) 次多项式（函数）。像这种用点表示多项式的手段我们叫做点值表达
噫，幼稚
我们先假设 \(f,g\) 次数相同都为 \(n\) ，那么我们的加法卷积 \(c\) 的次数为 \(2n\) ，换成点值表达就需要 \(2n+1\) 个点
于是幼稚的人类决定把系数表达换成点值表达，愚蠢的他们把这种操作成为 \(\mathtt{DFT}\)
但是他们并不喜欢用点值表达多项式，于是他们又想到了 \(\mathtt{DFT}\) 的逆运算，即将点值表达换成系数表达，叫做 \(\mathtt{IDFT}\)
乍一看麻烦了许多？
3. 复数与单位根
回顾一下 \(\mathtt{DFT}\) 的操作
找一堆数，代入进多项式得到点值
显然这些数是可以由我们自己决定的
那么有没有集好算与优秀与一身的数呢？
当然是有的
愚蠢的人类们的直觉是代入一些正整数啊，\(0\) 啊什么的进去，但是这些数并没有什么优秀的性质
直到一个聪明一点的人类跑出来，他叫傅里叶，是研究复数的专家（不知道复数的回炉重造去）
他尝试把复单位根代入进了多项式，发掘了一些伟大的性质
于是虫子们就掌握了 \(\mathtt{FFT}\)
所以，单位根是啥？
首先要从小学二年级数学老师跟你们讲的复数相乘讲起，还记得复数相乘在复平面的几何意义吗？
没错，模长相乘，幅角（逆时针）相加
\(n\) 次单位根是方程 \(x^n=1\) 的复数解
这俩有啥关系？
我们在复平面上画一个单位圆（借用某巨佬的图片）：
我们发现 \(1\) 这个数也在圆上（在复平面的圆是一堆数）
这些数的模长都为 \(1\)
换句话说，\(n\) 次单位根的解其实就在这个圆里，我们只需要在这个圆里面找单位根即可
显然幅角为 \(\frac{i}{n},i\in\{0,1,...,n-1\}\) 的数为 \(n\) 次单位根
为啥？考虑 \(x^n\) 的集合意义，即 \(x\) 这个数的模长到 \(|x|^n\) ，并且幅角为 \(n\times \text{x的幅角}\)
上述单位根经过旋转都到了 \(1\) 的位置
由于 \(n\) 次方程在复数域内有且仅有 \(n\) 个根，所以我们就不重不漏地找到了 \(n\) 次单位根
虫子数学家们把单位根叫做 \(\omega_{n}^{i},\omega_{n}^{0}=1\)， \(\omega_{n}^{i}\) 为第 \(i\) 个单位根
好奇的人类发现了关于单位根的很多性质：
- \(\omega_{n}^{k}=(\omega_{n}^{1})^k\)
- \(\omega_{n}^{i}\times\omega_{n}^{j}=\omega_{n}^{i+j}\)
- \(\omega_{2n}^{2k}=\omega_{n}^{k}\)
- \(\omega_{n}^{k+n/2}=-\omega_{n}^{k}\)
上面的性质都非常简单，几何性质考虑即可
4. \(\mathtt{DFT}\) 的实现过程
还记得发现复数根性质的虫子数学家傅里叶吗？
知道了这个性质，人类暴力地把多项式补成 \(2\) 的幂次次数的形式，再把它分成奇数次和偶数次
\[f(x)=\sum\limits_{k=0}^{n}f[k]x^k \]\[fl(x)=\sum\limits_{k=0}^{n/2-1}F[2k]x^k \]\[fr(x)=\sum\limits_{k=0}^{n/2-1}F[2k+1]x^k \]于是显然有 \(f(x)=fl(x^2)+xfr(x^2)\)
尝试把单位根代入
\[f(\omega_{n}^{k})=fl((\omega_{n}^{k})^2)+\omega_{n}^{k}fr((\omega_{n}^{k})^2) \]用上面的性质进行一系列化简会发现：
\[\begin{cases}f(\omega_{n}^{k})=fl(\omega_{n/2}^{k})+\omega_{n}^{k}fr(\omega_{n/2}^{k})\\f(\omega_{n}^{k+n/2})=fl(\omega_{n/2}^{k})-\omega_{n}^{k}fr(\omega_{n/2}^{k})\end{cases} \]然后就一个正负号的区别，虫子们通过分治可以获得 \(f(x)\) 的点值表示
虫子们又在思考：如何求出单位根呢？
根据前面的性质 \(\omega_{n}^{k}=(\omega_{n}^{1})^k\) ，所以只需要求出 \(\omega_{n}^{1}\) 即可
如上图，\(\alpha\) 为 \(\frac{2\pi}{n}\)，显然有 \(B(\cos(\frac{2\pi}{n}),\sin(\frac{2\pi}{n}))\) ，这个就是 \(\omega_{n}^{1}\) 的实部与虚部（弧度表示）
然后把 \(\omega_{n}^{1}\) 乘 \(n\) 次得到 \(\{\omega_{n}^{0},\omega_{n}^{1},...,\omega_{n}^{n-1}\}\)
虫子：简单的昂~
5. \(\mathtt{IDFT}\) 的实现过程
先给出虫子的结论：把 \(\mathtt{DFT}\) 的 \(\omega_{n}^{1}\) 换成 \(\omega_{n}^{-1}\) ，把乘换成除即可
为啥？单位根反演！
令点值数列 \(G=DFT(F)\)
\[G[k]=\sum\limits_{i=0}^{n-1}(\omega_{n}^{k})^iF[i] \]直接繁衍反演：
然后 \(\omega_{n}^{-1}=\omega_{n}^{n-1}=\cos(\frac{2\pi}{n})-\sin(\frac{2\pi}{n})i\) ，代入 \(\omega_{0}^{n},\omega_{-1}^{n},...,\omega_{0}^{-n+1}\)
6. show me your code
实现出来虫子们发现 \(\mathtt{DFT,IDFT}\) 只有一个符号的区别（想必聪明的你已经发现了）
于是可以凑在一起
const double pi = acos(-1);
const int maxn = 2e6 + 200;
struct complexs {// 小学二年级小孩子都能用脚指头上的指甲盖上的蛋白质上的氨基酸想出来的复数结构体
double x, y;
complexs(double xx = 0, double yy = 0) { x = xx, y = yy; }
complexs operator + (complexs const &rhs) const { return complexs(x + rhs.x, y + rhs.y); }// 复数相加
complexs operator - (complexs const &rhs) const { return complexs(x - rhs.x, y - rhs.y); }// 复数相减
complexs operator * (complexs const &rhs) const { return complexs(x * rhs.x - y * rhs.y, x * rhs.y + y * rhs.x); }// 复数相乘
complexs operator / (complexs const &rhs) const { return complexs((x * rhs.x + y * rhs.y) / (rhs.x * rhs.x + rhs.y * rhs.y), (y * rhs.x - x * rhs.y) / (rhs.x * rhs.x + rhs.y * rhs.y)); }// 复数相除
} f[maxn << 1], p[maxn << 1], sav[maxn << 1];
int n, m;
void fft(complexs *f, int len, bool flag) {
if (len == 1) return;// 如果分治到头了，虫子返回
complexs *fl = f, *fr = f + len / 2;// 得到两个系数/点值数列
for (int i = 0; i < len; i++) sav[i] = f[i];// copy一份数列
for (int i = 0; i < len / 2; i++) {// 处理奇偶次数的数列
fl[i] = sav[i << 1];
fr[i] = sav[i << 1 | 1];
}
// 虫子都会的分治
fft(fl, len / 2, flag);
fft(fr, len / 2, flag);
// 得到前面说的第二个单位根和第一个单位根
complexs base(cos(2 * pi / len), sin(2 * pi / len)), p(1, 0);
base.y *= (flag ? 1 : -1);// DFT or IDFT
for (int i = 0; i < len / 2; i++) {
sav[i] = fl[i] + p * fr[i];// 根据上面的柿子计算左边
sav[i + len / 2] = fl[i] - p * fr[i];//根据上面的柿子计算右边
p = p * base;// 下一个单位根
}
for (int i = 0; i < len; i++) f[i] = sav[i];// copy回来
}
int main() {
n = read();
m = read();
for (int i = 0; i <= n; i++) f[i].x = read();
for (int i = 0; i <= m; i++) p[i].x = read();
for (m += n, n = 1; n <= m; n <<= 1) {}// 补全次数
fft(f, n, 1);// 分治 DFT 处理 f
fft(p, n, 1);// 分治 DFT 处理 p
for (int i = 0; i < n; i++) f[i] = f[i] * p[i];// 将两个点值表达乘起来
fft(f, n, 0);// 搞回系数表达
for (int i = 0; i <= m; i++) {
write((int)(f[i].x / n + 0.49));// 单位根反演的除以 n ，四舍五入
putchar(' ');
}
return 0;
}
\(\mathtt {NTT}\) 部分
概述
由于 \(\mathtt{FFT}\) 使用了三角函数之类的计算，所以精度未免误差太大
虫子们又在想，有没有什么办法能够完美地避开这个误差呢？
想想，虫子们是怎么精确计算分数的？
没错，取模、乘法逆元
所以我们需要在模意义下找到那么一个东西，满足单位根的性质
很明显这东西和原根有关：
众虫周知，对于素数 \(p\) ，其模意义下的同余系恰恰有 \(1,2,3,...,p-1\) 这么多数
对于 \(g\) ，称 \(g\) 为 \(p\) 的原根，当且仅当 \(g\) 的阶达到 \(p-1\) 的上界
显然 \(p-1\) 并不能等于 \(n\) ，所以原根本身还不能用作单位根
但是聪明的虫子们能够发现 \(g^k\) 的阶数为 \(\frac{p-1}{\gcd(p-1, k)}\) ，这个数为 \(p-1\) 的约数
当 \(n\!\not|\ (p-1)\) 时，找不到阶恰为 \(n\) 的数
当 \(n\mid (p-1)\) 时，不难发现 \(g^{\frac{p-1}{n}}\) 的阶恰为 \(n\)
这说明 \(g^{\frac{p-1}{n}}\) 满足了单位根的前 \(3\) 个条件，第 \(4\) 个也易证
由于满足 \(n\mid(p-1)\) 才能构造出单位根，而在上文中虫子们把 \(n\) 补成了 \(2\) 的幂，所以现在只需要找出一个含有大量 \(2\) 的幂的数 \(p-1\) 即可
有一个熟悉的数满足这个条件：
\(998244353=2^{23}\times 7\times 17+1\)
所以虫子们一般用 \(p=998244353\) 来计算 \(\mathtt{NTT}\)
代码略有变化，建议全文背诵
const int maxn = 3e6 + 300;
const int mod = 998244353;
const int gi = 332748118;
const int g = 3;
int n, m, inv, len, lim = 1, tr[maxn];
int f[maxn], p[maxn];
int ksm(int p, int q) {
int res = 1;
while (q) {
if (q & 1) res = res * p % mod;
p = p * p % mod;
q >>= 1;
}
return res;
}
inline void ntt(int *f, int flag){
for (int i = 0; i < lim; i++) {
if (i < tr[i]) swap(f[i], f[tr[i]]);
}
for (int p = 1; p < lim; p <<= 1) {
int tg = ksm(flag == 1 ? g : gi, (mod - 1) / (p << 1));
for (int j = 0, w = 1; j < lim; j += (p << 1), w = 1) {
for (int k = 0; k < p; k++, w = (w * tg) % mod) {
int x = f[j + k];
int y = w * f[j + k + p] % mod;
f[j + k] = (x + y) % mod;
f[j + k + p] = (x - y + mod) % mod;
}
}
}
}
signed main() {
n = read();
m = read();
for (int i = 0; i <= n; i++) f[i] = (read() + mod) % mod;
for (int i = 0; i <= m; i++) p[i] = (read() + mod) % mod;
while (lim <= n + m) lim <<= 1, len++;
inv = ksm(lim, mod - 2);
for (int i = 0; i < lim; i++) tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << (len - 1));
ntt(f, 1);
ntt(p, 1);
for (int i = 0; i < lim; i++) f[i] = (f[i] * p[i]) % mod;
ntt(f, -1);
for (int i = 0; i <= n + m; i++) {
write((f[i] * inv) % mod);
putchar(' ');
}
return 0;
}
标签：const,int,mathtt,rhs,FFT,单位根,NTT,笔记,omega 来源： https://www.cnblogs.com/Ender32k/p/15780957.html