FFT / NTT 入门笔记
前言
感觉如果认真学真的很快就能学懂。（？）
给自己看的，可能比较简略，如果有不懂的可以去看看后记的参考资料。
其中全家桶的顺序是按照我自己学的顺序写的，其实可能 ln 放在最前面可能更合适。
引入
考虑以下问题：
给定 \(f_{0\sim n}, g_{0\sim m}\)。
定义 \(f\times g = h\) 表示 \(\forall i\in [0, n + m], h_i = \sum\limits_{j = \max(0, i - m)}^{\min(n, i)} f_j g_{i - j}\)。
现在需要给出 \(h_{0\sim n + m}\) 的值。
如果直接按照定义式求解，时间复杂度将会是 \(\mathcal{O}(n^2)\) 的。
而 FFT, NTT 这两个算法可以将求解该问题的复杂度优化到 \(\mathcal{O}(n\log n)\)。
FFT & NTT
多项式与点值
首先换一下这个 \(f\times g = h\) 的形式：\(\forall i\in [0, n], j\in [0, m], h_{i + j}\gets f_i g_j\)。
考虑这个 \(i, j \to i + j\) 的形式并不契合 \(f_i, g_j \to f_i g_j\) 的形式，尝试用一个占位符去同时表示这个 \(f_i\times g_j\) 和 \(i + j\) 的变化。
对此考虑到 \(x^i\times x^j = x^{i + j}\)，考虑把 \(f, g, h\) 都转化为多项式的形式：\(F(x) = \sum\limits_{i = 0}^n f_i x^i\)，对于 \(G(x), H(x)\) 也同理，同时记 \([x^i]F(x)\) 就是占位符 \(x^i\) 的系数，也就对应 \(f_i\)。
那么此时 \(f\times g = h\) 就直接对应多项式乘法 \(F(x)\times G(x) = H(x)\)，也称作卷积。
令多项式 \(F(x)\) 的最高次为 \(f_i\) 不为 \(0\) 的 \(i\) 的最大值。
接下来介绍点值：
点值其实就是把 \(F(x)\) 中的 \(x\) 代入一个数值 \(x_i\) 带进这个求和得到的值，记 \(y_i = F(x_i)\)。
记 \(F(x)\) 的最高次为 \(n\)，那么只需要知道 \(n + 1\) 个不同的点值 \((x_i, y_i)(\forall 1\le i < j\le n + 1, x_i\not = x_j)\)，就可以唯一确定这个多项式的 \(n + 1\) 个系数（\(f_{0\sim n}\)）。
证明：
考虑把这个代入点值的形式写作矩阵，有 \(\begin{bmatrix}f_0& f_1& \cdots & f_n \end{bmatrix}\times \begin{bmatrix}1 &x_1 & x_1^2 & \cdots & x_1^n\\ 1 & x_2 & x_2^2 & \cdots & x_2^n\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n + 1} & x_{n + 1}^2 & \cdots & x_{n + 1}^n\end{bmatrix} = \begin{bmatrix}y_1& y_2& \cdots & y_{n + 1}\end{bmatrix}\)，记做 \(F\times X = Y\)。
那么只要 \(X\) 存在逆 \(X^{-1}\)，就有 \(Y\times X^{-1} = F\)，就能说明能够通过 \(Y\) 推出 \(F\) 了。
而 \(X\) 是范德蒙德矩阵，有 \(\det(X) = \prod\limits_{1 \le i < j\le n + 1} (x_j - x_i)\not = 0\)，所以存在 \(X^{-1}\)。
并且考虑到如果 \(F(x) \times G(x) = H(x)\)，那么 \(F(x_i)\times G(x_i) = H(x_i)\)。
于是假设 \(H(x)\) 的最高次为 \(n\)，一个想法是：
- 在 \(F(x)\) 和 \(G(x)\) 代入 \(n + 1\) 个不同的 \(x_i\) 得到点值 \((x_i, F(x_i)), (x_i, G(x_i))\)。
- 因为 \(F(x_i)\times G(x_i) = H(x_i)\)，只需对应相乘便可得到 \(n + 1\) 个不同的点值 \((x_i, H(x_i))\)。
- 用 \(n + 1\) 个不同的点值 \((x_i, H(x_i))\) 去确定 \(H(x)\) 的系数。
其中第二步的复杂度是 \(\mathcal{O}(n)\) 的，不再考虑。
称第一步为 DFT，第三步为 IDFT。
单位根（原根）与 DFT
引入 FFT-DFT 所需要的单位根这个概念：
考虑 \(x^n = 1\) 这个方程，这个方程会存在 \(n\) 个复数意义下的根。
这 \(n\) 个根分别是 \(\omega_n^k = \cos(\frac{2\pi}{n}\times k) + \sin(\frac{2\pi}{n}\times k)i(0 \le k < n)\)（在实轴及虚轴所构成的二维平面内把单位圆 \(n\) 等分，每一部分对应的值），这里的 \(i\) 指的是虚数单位。
而 \(\omega_n^k\) 有着一些特别好的性质：
- \(\omega_n^i\times \omega_n^j = \omega_n^{i + j}\)。
- \(\omega_n^k = \omega_n^{k + n} = \omega_n^{k - n}\)。（但是本质不同的还是只有 \(n\) 个，所以放到 \([0, n)\) 来考虑）。
这说明可以把 \(k\) 放在模 \(n\) 意义下来看。 - 若 \(n\bmod 2 = 0\)，则 \(\omega_n^k = -\omega_n^{k + \frac{n}{2}}\)。
- 若 \(n\bmod 2 = 0, k\bmod 2 = 0\)，则 \(\omega_n^k = \omega_{\frac{n}{2}}^{\frac{k}{2}}\)。
考虑把 \(H(x)\) 的最高次 \(n\)（\(F(x), G(x)\) 的最高次之和）手动补足到 \(2^m - 1\)， 并把 \(F(x), G(x), H(x)\) 填充 \(0\) 至最高次也为 \(2^m - 1\)。
那么令 \(n + 1\) 个点值就为 \((\omega_n^k, F(\omega_n^k))(0\le k < n)\)，考虑上述形式的对于这个求解的帮助。
考虑要求得 \(\sum\limits_{i = 0}^{n - 1} f_i \omega_n^{ki}\)，那么展开来写：\(f_0 \omega_n^0 + f_1\omega_n^k + \cdots + f_{n - 2}\omega_n^{k(n - 2)} + f_{n - 1}\omega_n^{k(n - 1)}\)。
考虑让形式尽量一致，于是分开奇偶项：\((f_0\omega_n^0 + \cdots + f_{n - 2}\omega_n^{k(n - 2)}) + \omega_n^k(f_1 \omega_n^0 + \cdots + f_{n - 1}\omega_n^{k(n - 2)})\)。
此时括号内的指标都应该为偶数，于是考虑再化简：\((f_0\omega_{\frac{n}{2}}^0 + \cdots + f_{n - 2}\omega_{\frac{n}{2}}^{\frac{k(n - 2)}{2}}) + \omega_n^k(f_1\omega_{\frac{n}{2}}^0 + \cdots + f_{n - 1}\omega_{\frac{n}{2}}^{\frac{k(n - 2)}{2}})\)。
记 \(FL(x) = f_0x^0 + \cdots + f_{n - 2}x^{\frac{n - 2}{2}}, FR(x) = f_1x^0 + \cdots + f_{n - 1}x^{\frac{n - 2}{2}}\)。
于是会发现如果知道了 \(FL(\omega_{\frac{n}{2}}^k), FR(\omega_{\frac{n}{2}}^k)\)，就可以得到 \(F(\omega_n^k) = FL(\omega_\frac{n}{2}^k) + \omega_n^kFR(\omega_{\frac{n}{2}}^k)\)。
看起来好像对于每一个 \(k\) 都需要知道 \(FL(\omega_\frac{n}{2}^k), FR(\omega_\frac{n}{2}^k)\) 的值，但是再考虑：\(\omega_\frac{n}{2}^k = \omega_\frac{n}{2}^{k + \frac{n}{2}}\)。
这说明只需要知道 \(0\le k < \frac{n}{2}\) 中 \(FL(\omega_\frac{n}{2}^k), FR(\omega_\frac{n}{2}^k)\) 的值即可。
那么此时对于长度为 \(\frac{n}{2}\) 的 \(FL(x), FR(x)\) 也只需要算 \(\forall 0\le k < \frac{n}{2}, \omega_\frac{n}{2}^k\) 的点值，且合并出 \(F(\omega_n^k)\) 是 \(\mathcal{O}(n)\) 的（每个值都可以直接求）。
所以总的时间复杂度是 \(T(n) = 2T(\frac{n}{2}) + \mathcal{O}(n)\)，有 \(T(n) = \mathcal{O}(n\log n)\)。
一些减少常数的技巧：
- 发现对于 \(0\le k < \frac{n}{2}\)，\(F(\omega_n^k) = FL(\omega_\frac{n}{2}^k) + \omega_n^kFR(\omega_{\frac{n}{2}}^k), F(\omega_n^{k + \frac{n}{2}}) = FL(\omega_\frac{n}{2}^k) - \omega_n^kFR(\omega_{\frac{n}{2}}^k)\)。
所以可以直接用 \(FL(\omega_\frac{n}{2}^k), \omega_n^k FR(\omega_\frac{n}{2}^k)\) 这两个值一次性求出 \(F(\omega_n^k), F(\omega_n^{k + \frac{n}{2}})\)。 - 把递归形式优化成迭代的形式。
那么主要问题是把 \(F\) 分离成 \(FL, FR\) 的问题。
观察会发现这个过程就是类似于，把偶数的位放到前面，奇数的位放在后面，然后两堆分别抹掉最低位继续递归下去。
于是对于 \(i\) 这个位置在递归到叶子的位置时在的位置应该是 \(i\) 二进制翻转后的值，可以提前预处理。
这也被称为“蝴蝶变换”。
对于 NTT-DFT 其实更为常用，因为这类问题通常都是用来解决计数问题。
而这类问题通常都需要模数，而其中的一些模数因为其性质就适合 NTT。
接下来引入原根：
因为对于质数 \(p\)，\(\varphi(p) = p - 1\)，关于质数 \(p\) 的原根 \(g\) 需要满足 \(g^{p - 1}\equiv 1\pmod p\)，\(\forall i\in [1, p - 1), g^i \not \equiv 1\pmod p\)。
十分常用的就是 \(p = 998244353, g = 3\)，这也是因为其独特的性质：
- \(p = 998244353 = 119\times 2^{23} + 1\)，于是 \(\varphi(p) = 119\times 2^{23}\)。
- 对于 \(g^{119\times 2^i}(0\le i\le 23)\)，其有着模 \(p\) 意义下 \(2^{23 - i}\) 的循环节。
也即若想要有模 \(p\) 意义下 \(2^i\) 的循环节，可以选取 \(g^{119\times 2^{23 - i}} = g^{\frac{\varphi(p)}{2^i}}\)。 - \(g^{\frac{\varphi(p)}{2^i}}(0\le i\le 23)\) 在模 \(p\) 意义下满足前面单位根所需要的所有性质。（\(\omega_{2^i}^k\) 可以看作 \(g^{\frac{\varphi(p)}{2^i}\times k}\)）。
也就是说 \(p = 998244353\) 可以解决总长度在 \(2^{23}\) 以内的多项式乘问题。
NTT-DFT 一些减小常数的技巧：
- 每迭代一次就会在每个值上加一个 \(\le p\) 的值，如果使用
u64可以迭代到处理完 \(2^{17}\) 时再对所有 \(F(x_i)\) 取模一次（\(18\times 998244353^2 < 2^{64}\)）。
不过这个指的是求得 \(F(x)\) 的时候涉及的加减法不用取模，中间乘法的步骤还是需要取模的。 - FFT-DFT 的优化也可以用到这里。
值得一提的是，因为 FFT 需要维护实部虚部都为浮点数的复数，NTT 的实际效率会比较优秀。
IDFT
其实说起来十分简单，对于 FFT-IDFT，步骤与 FFT-DFT 基本相同，只需要把单位元 \(\omega_{n}^1\) 替换成 \(\omega_n^{-1}\) 求解，最后每一项的系数都除掉项数 \(n\)（注意这里不是最高次了）。
对于 NTT-IDFT 也是类似的，只需要把 NTT-DFT 中的 \(g\) 换乘 \(g^{-1}\) 求解，最后系数除掉项数 \(n\) 即可。
证明考虑回到说明 \(n + 1\) 个不同点值可以确定 \(n\) 次多项式那里。
相当于是证明：
\[\begin{bmatrix}\omega_n^{0\times 0} & \omega_n^{0\times 1} & \cdots & \omega_n^{0\times (n - 1)}\\ \omega_n^{1\times 0} & \omega_n^{1\times 1} & \cdots & \omega_n^{1\times (n - 1)}\\ \vdots & \vdots & \ddots & \vdots \\ \omega_n^{(n - 1)\times 0} & \omega_n^{(n - 1)\times 1} & \cdots & \omega_n^{(n - 1)\times (n - 1)}\end{bmatrix}\times \begin{bmatrix}\omega_n^{-0\times 0} & \omega_n^{-0\times 1} & \cdots & \omega_n^{-0\times (n - 1)}\\ \omega_n^{-1\times 0} & \omega_n^{-1\times 1} & \cdots & \omega_n^{-1\times (n - 1)}\\ \vdots & \vdots & \ddots & \vdots \\ \omega_n^{-(n - 1)\times 0} & \omega_n^{-(n - 1)\times 1} & \cdots & \omega_n^{-(n - 1)\times (n - 1)}\end{bmatrix} = \begin{bmatrix}n & 0 & \cdots & 0\\ 0 & n & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & n\end{bmatrix} \]即 DFT 矩阵和 IDFT 矩阵在相乘后为单位矩阵 \(\times n\)，也就是 DFT 矩阵与 IDFT 矩阵 \(\times \frac{1}{n}\) 互为逆矩阵。
考虑直接代入矩阵乘法的公式：\(v_{i, j} = \sum\limits_{k = 0}^{n - 1} \omega_n^{ik}\omega_n^{-kj} = \sum\limits_{k = 0}^{n - 1} \omega_n^{k(i - j)}\)，然后分讨：
- \(i = j\)，那么 \(v_{i, j} = \sum\limits_{k = 0}^{n - 1} \omega_n^0 = n\)。
- \(i\not = j\)，那么 \(v_{i, j} = \frac{\omega_n^{n(i - j)} - 1}{\omega_n^{i - j} - 1} = \frac{1 - 1}{\omega_n^{i - j} - 1} = 0\)。
多项式杂项 & 全家桶
多项式加减法与乘常数
加减法就是对于相同的占位符 \(x^i\) 的系数相加，即 \(F(x) \pm G(x) = H(x)\) 就有 \([x^i]H(x) = [x^i]F(x)\pm G(x)(i\ge 0)\)。
乘常数就是对于每一个占位符的系数都乘上该常数，即 \(cF(x) = G(x)\) 就有 \([x^i]G(x) = c\times [x^i]F(x)\)。
时间复杂度均为 \(\mathcal{O}(n)\)。
多项式求导与积分
令 \(F'(x) = G(x)\)，那么有 \(G(x) = \displaystyle\sum\limits_{i = 1} i\times [x^i]F(x)\times x^{i - 1}\)，\(F(x) = C + \displaystyle\sum\limits_{i = 0} \frac{1}{i + 1}\times [x^{i}]G(x) \times x^{i + 1}\)。
其中 \(F(x)\to G(x)\) 就是求导的过程，\(G(x)\to F(x)\) 就是积分的过程。
时间复杂度均为 \(\mathcal{O}(n)\)。
多个多项式卷积
若现在需要得到 \(A(x) \times B(x)\times C(x) = D(x)\)。
一个想法是先卷 \(A(x)\times B(x) = E(x)\)，再卷 \(E(x) \times C(x) = D(x)\)。
不过这样就需要 \(3\) 次 DFT，\(2\) 次 IDFT 了。
考虑到因为点值是可以直接相乘的，就算现在是 \((x_i, A(x_i)), (x_i, B(x_i)), (x_i, C(x_i))\)，也照样有 \(A(x_i)B(x_i)C(x_i) = D(x_i)\)。
所以在长度允许的情况下可以直接 \(3\) 次 DFT 得到 \(A(x), B(x), C(x)\) 的点值，再一次性相乘得到 \(D(x)\) 的点值，这样就只需要 \(1\) 次 IDFT 得到 \(D(x)\) 了。
类似的，当有 \(m\) 个多项式卷积时，在长度允许的情况下也可以只用 \(m\) 次 DFT 和 \(1\) 次 IDFT 解决。
如果长度并不在允许的范围内，可以尝试分组合并，即通过分组让每一组的卷积长度在允许范围内，再对每一组得到的多项式继续分组做下去直到只剩 \(1\) 个。
三次变两次优化
考虑 FFT 的时候虽然扩宽到了复数域，也就多了个虚部要维护。
但是在卷积过程中，初始的 \(F(x), G(x)\) 和结果 \(H(x)\) 的虚部都为 \(0\)，那能不能把虚部也利用上呢。
于是尝试考虑在虚部处构造出 \(F(x)G(x)\)，也即让最后卷出来的结果带 \(F(x)G(x)i\)。
此时考虑完全平方公式：\((a + b)^2 = a^2 + b^2 + 2ab\)。如果把 \(b\) 代入为 \(bi\) 就有：\((a + bi)^2 = a^2 - b^2 + 2abi\)。
发现这个 \(2abi\) 就完全满足了想要的 \(F(x)G(x)i\) 这一个样子。
于是考虑令 \(H(x) = F(x) + G(x)i\)，考察 \(H^2(x) = (F^2(x) - G^2(x)) + 2F(x)G(x)i\)，最后虚部 \(\times \frac{1}{2}\) 后的多项式就是想要的结果。
可以考虑类似上问的多个多项式卷积，得到 \((x_i, H(x_i))\) 后，\((H(x_i))^2 = H^2(x_i)\)，也就直接得到了 \(H^2(x)\) 的点值 \((x_i, H^2(x_i))\)。
于是这只需要 \(1\) 次 FFT-DFT 和 \(1\) 次 FFT-IDFT，省掉了 \(1\) 次 DFT。
至于在更多的多项式卷积时还能否有妙用，我并没有探究出来。
看起来这个想法其实是把 \(i^k\) 也当作了占位符，但是 \(k\) 是在模 \(2\) 意义下的，当有 \(\ge 3\) 个多项式时，难免会与 \(i^{k - 2}, i^{k + 2}\) 等产生一些冲突，而且只能放置在 \(i^0, i^1\) 下也限制了一定的形式。
例如若有 \(3\) 个时构造成 \(F(x) + (G(x) + H(x))i\)，在 \(3\) 次方下 \(i^0, i^2\) 所对应的系数是 \(F^3(x) - 6F(x)G(x)H(x) - 3F(x)G^2(x) - 3F(x)H^2(x)\)，干扰的系数看起来并不好处理。
牛顿迭代处理 inv, sqrt, exp
首先需要用到泰勒展开：
对于 \(F(x)\) 考察在 \(a\) 处的展开有：\(F(x) = \sum\limits_{i = 0}\dfrac{F^{(i)}(a)}{i!}(x - a)^i\)，大致就是一个以 \(a\) 这里展开并无限逼近的过程。
牛顿迭代处理多项式的一些变换的步骤是：
- 构造出 \(G(F(x))\) 要求得满足 \(G(F(x)) \equiv 0\pmod {x^n}\) 的 \(F(x)\bmod x^n\)。
- 直接算出常数项的值。
- 不妨设 \(n\) 是偶数，通过满足 \(G(F_0(x)) \equiv 0\pmod {x^{\frac{n}{2}}}\) 的 \(F_0(x)\bmod x^\frac{n}{2}\) 来进行操作得到 \(F(x)\)。
第一步第二步的构造待具体例子再提，先就考虑第三步。
首先 \(G(F(x)) \equiv 0\pmod {x^n}\) 那肯定也就有 \(G(F(x))\equiv 0\pmod {x^\frac{n}{2}}\)。
再考虑 \(G(F(x)) = \sum\limits_{i = 0} [x^i]G(x)\times F^i(x)\)，那么在 \(\bmod x^\frac{n}{2}\) 意义下就只会有 \(F(x) \bmod x^\frac{n}{2}\) 有用。
又因为 \(F_0(x)\bmod x^\frac{n}{2}\) 就已经满足了 \(G(F_0(x))\equiv 0\pmod {x^\frac{n}{2}}\)，所以可以直接令 \([x^i]F(x) = [x^i]F_0(x)(0\le i < \frac{n}{2})\)，即直接让 \(F(x)\) 的低 \(\frac{n}{2}\) 位就为 \(F_0(x)\)。
此时再考察 \(G(F(x))\) 在 \(F_0(x)\) 处的展开有：\(G(F(x)) = \sum\limits_{i = 0}\dfrac{G^{(i)}(F_0(x))}{i!}\left(F(x) - F_0(x)\right)^i\)。
又因为 \(G(F(x))\equiv 0\pmod {x^n}\)，所以有 \(\sum\limits_{i = 0}\dfrac{G^{(i)}(F_0(x))}{i!}\left(F(x) - F_0(x)\right)^i \equiv 0 \pmod {x^n}\)。
再考虑到因为 \(F(x)\) 的低 \(\frac{n}{2}\) 位与 \(F_0(x)\) 相同，所以 \(F(x) - F_0(x)\equiv 0\pmod {x^\frac{n}{2}}\)，那么就有 \((F(x) - F_0(x))^\alpha \equiv 0\pmod {x^n}(\alpha \ge 2)\)。
所以上述的展开式只在 \(i = 0, 1\) 有用，即 \(G(F_0(x)) + G'(F_0(x))(F(x) - F_0(x)) \equiv 0\pmod {x^n}\)。
于是可以解得 \(F(x) \equiv F_0(x) - \dfrac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n}\)，于是便可以通过多项式的其他变换进行求解 \(F(x)\)。
一些可能的常数优化：
- 把递归求 \(F_0(x)\) 换为倍增。
接下来就是多项式 inv, sqrt, exp 在牛顿迭代下的推导：
-
多项式 inv，即给定 \(H(x)\)，要求 \(H^{-1}(x) \equiv F(x)\pmod {x^n}\)。
那么可以构造 \(G(F(x)) = F^{-1}(x) - H(x)\)，对于常数项就是 \([x^0]F(x) = ([x^0] H(x))^{-1}\)，需要保证 \([x^0]H(x)\not = 0\)。
代入 \(G(F(x))\)，有 \(F(x) \equiv F_0(x) - \dfrac{F_0^{-1}(x) - H(x)}{-F_0^{-2}(x)} \equiv 2F_0(x) - F_0^2(x)H(x) \equiv F_0(x)(2 - F_0(x)H(x))\pmod {x^n}\)。
那么只需要用多项式乘法就可以从 \(F_0(x)\) 得到 \(F(x)\)。时间复杂度为 \(T(n) = T(\frac{n}{2}) + \mathcal{O}(n\log n)\)，那么 \(T(n) = \mathcal{O}(n\log n)\)。
-
多项式 sqrt，即给定 \(H(x)\)，要求 \(H^{\frac{1}{2}}(x) \equiv F(x)\pmod {x^n}\)。
那么可以构造 \(G(F(x)) = F^2(x) - H(x)\)，对于常数项就是 \([x^0]F(x) = ([x^0]H(x))^{\frac{1}{2}}\)，需要求一个二次剩余，但是如果 \([x^0]H(x) = 1\) 就可以直接用 \(1\) 啦。
代入 \(G(F(x))\)，有 \(F(x) \equiv F_0(x) - \dfrac{F_0^2(x) - H(x)}{2F_0(x)} \equiv \dfrac{F_0^2(x) + H(x)}{2F_0(x)}\pmod {x^n}\)。
那么用多项式乘法和 inv 就可以从 \(F_0(x)\) 得到 \(F(x)\)。时间复杂度分析同上为 \(\mathcal{O}(n\log n)\)。
-
多项式 exp，即给定 \(H(x)\)，要求 \(\exp(H(x)) \equiv F(x)\pmod {x^n}\)。
那么可以构造 \(G(F(x)) = \ln F(x) - H(x)\)，对于常数项就要求 \([x^0]H(x) = 0\)，有 \([x^0]F(x) = 1\)。
代入 \(G(F(x))\)，有 \(F(x)\equiv F_0(x) - \dfrac{\ln F_0(x) - H(x)}{\frac{1}{F_0(x)}}\equiv F_0(x)(1 - \ln F_0(x) + H(x))\pmod {x^n}\)。
那么用多项式乘法和 ln 就可以从 \(F_0(x)\) 推到 \(F(x)\)，不过现在还没有说到 ln，所以暂时还写不了。
带余除法
给定 \(F(x), G(x)\)，记 \(n\) 为 \(F(x)\) 的项数，\(m\) 为 \(G(x)\) 的项数，要求 \(Q(x), R(x)\) 满足 \(F(x) = G(x)Q(x) + R(x)\)，且 \(R(x)\) 的项数 \(< m\)（如果 \(\ge\) 就可以再除下去）。
不妨直接令 \(R(x)\) 项数为 \(m - 1\)，那么只会填充一些 \(0\) 并无影响。
一个非常妙的想法是，考虑如果可以先舍弃掉 \(R(x)\) 那么就是一个多项式 inv 和乘法，所以尝试去抛掉 \(R(x)\)。
记一个多项式 \(A(x)\) 最高次为 \(n_A\)，定义 \(A_R(x) = x^{n_A}A(x^{-1})\)，其实也就是把 \(A\) 的系数翻转了。
考虑 \(F(x) = G(x)Q(x) + R(x)\)，那么就有 \(F(x^{-1}) = G(x^{-1})Q(x^{-1}) + R(x^{-1})\)。
以下需要注意最高次等于项数 \(-1\)，于是有 \(x^{n - 1}F(x^{-1}) = x^{n - 1}G(x^{-1})Q(x^{-1}) + x^{n - 1}R(x^{-1})\)，那么换成 \(F_R(x) = G_R(x)Q_R(x) + x^{n - m + 1}R_R(x)\)。
回到最初的想法：想要抛掉 \(R(x)\)。
于是考虑模 \(x^{n - m + 1}\)，那就有 \(F_R(x)\equiv G_R(x)Q_R(x)\pmod {x^{n - m + 1}}\)，那么就有 \(Q_R(x) \equiv \dfrac{F_R(x)}{G_R(x)}\pmod {x^{n - m + 1}}\)。
于是就可以通过多项式 inv 和乘法求出 \(Q_R(x)\)，那就得到了 \(Q(x)\)，此时再用多项式乘法得到 \(R(x) = F(x) - G(x)R(x)\)。
求 \(F_R(x)\) 因为就是翻转系数是 \(\mathcal{O}(n)\) 的，时间复杂度仍然是 \(\mathcal{O}(n\log n)\)。
ln, exp
给定 \(F(x)\)，要求 \(\ln F(x)\equiv G(x)\pmod {x^n}\)。
考虑直接在两侧求导，注意左边是个复合导，有 \(\dfrac{F'(x)}{F(x)} \equiv G(x)\pmod {x^n}\)。
于是只需要用多项式求导，inv，乘法和积分就可以了，时间复杂度 \(\mathcal{O}(n\log n)\)。
这里再提到 exp 的原因只是因为刚刚牛迭时还不会处理多项式 ln。
推导过程见上，这里直接给出最后结论，倍增时有 \(F(x)\equiv F_0(x)(1 - \ln F_0(x) + H(x))\pmod {x^n}\)。
那么用多项式乘法和 ln 就可以解决，时间复杂度分析也和其他牛迭处理的变换一样，为 \(\mathcal{O}(n\log n)\)。
pow
给定 \(F(x), k\)，要求 \(F^k(x) \equiv G(x)\pmod {x^n}\)。
考虑如果是求 \(a^b = c\)，一个办法是取 \(\ln\)，即 \(b\ln a = \ln c\) 得到 \(\ln c\)，那么再 \(\exp\) 回去就可以了。
考虑把同样的思想套用过来，那么就有 \(G(x) \equiv \exp(k\ln F(x))\pmod {x^n}\)。
不过这里还有一些疑问：
-
\(\ln F(x)\) 需要保证 \([x^0]F(x) = 1\)，那如果不满足怎么办呢？
首先对于 \([x^0]F(x) = 0\) 的情况，会发现此时 \([x^i]G(x)(0\le i < k)\) 一定都为 \(0\)（\(F(x)\) 最小有值位为 \(1\)，那么 \(G(x)\) 最小有值位就是 \(1\times k = k\)）。
所以此时直接舍弃 \(F(x)\) 的常数项，并舍弃 \(G(x)\) 前 \(k\) 项就行了。对于 \([x^0]F(x) > 1\) 的情况，那么让 \([x^0]F(x) = 1\) 不就行了吗。
所以记 \(v = [x^0]F(x)\)，那么就有 \(G(x) \equiv v^k \exp(k\ln (v^{-1}F(x)))\pmod {x^n}\)。 -
当 \(k\) 很大的时候，直接对 \(p\) 取模为什么是对的？
考虑泰勒展开：\(\exp (k\ln F(x)) = \displaystyle\sum\limits_{i = 0}\dfrac{(k\ln F(x))^i}{i!} = \displaystyle\sum\limits_{i = 0}\dfrac{k^i \ln^i F(x)}{i!}\)，便可以说明了。
实现时可能还需要注意一下 \(F(x)\) 系数全为 \(0\) 和 \(k = 0\) 的情况。
需要用到多项式 ln 和 exp，时间复杂度为 \(\mathcal{O}(n\log n)\)。
代码（NTT 全家桶）
注意事项：
- 开头需要调用
init()，其中LIM_N的值代表卷积的长度不会超过1 << LIM_N。 - 所有
n, m都代表多项式项数（或者说是最高次 +1）。 inv需要保证常数项不为 \(0\)。sqrt只实现了常数项为 \(1\) 的情况，二次剩余还没有写。ln需要保证常数项为 \(1\)。exp需要保证常数项为 \(0\)。pow需要传入 \(3\) 个参数：k1, k2, k3，分别代表 \(k\bmod p, k\bmod (p - 1), \min(k, L)(L > n)\)。
这三个分别用来处理 exp 时的 \(k\)，\(v^k\) 时的 \(k\)，以及移位时所需要的 \(k\)。
当 \(k\) 能直接用u64存下时也可以直接只传入一个 \(k\)。
#include <bits/stdc++.h>
using u64 = unsigned long long;
constexpr u64 mod = 998244353;
constexpr u64 G = 3, invG = 332748118;
constexpr u64 inv2 = (mod + 1) / 2;
constexpr inline u64 qpow(u64 a, u64 b) {
u64 v = 1;
for (; b; b >>= 1, a = a * a % mod) {
if (b & 1) v = v * a % mod;
}
return v;
}
inline void add(u64 &x, const u64 &y) { (x += y) >= mod && (x -= mod); }
namespace ntt {// ---------------------------------------------------------------------
using poly = std::vector<u64>;
constexpr int LIM_N = 18;
u64 rG[LIM_N + 1], riG[LIM_N + 1];
int rev[1 << LIM_N];
u64 invi[1 << LIM_N];
#define ctz __builtin_ctz
inline void init() {
for (int i = 1; i <= LIM_N; i++) rG[i] = qpow(G, (mod - 1) / (1 << i));
for (int i = 1; i <= LIM_N; i++) riG[i] = qpow(invG, (mod - 1) / (1 << i));
for (int i = 0; i < (1 << LIM_N); i++) {
rev[i] = rev[i >> 1] >> 1 | (i & 1) << LIM_N - 1;
}
invi[0] = invi[1] = 1;
for (int i = 2; i < (1 << LIM_N); i++) {
invi[i] = (mod - mod / i) * invi[mod % i] % mod;
}
}
inline poly dft(poly f, int n) {
const int shift = LIM_N - ctz(n);
for (int i = 0; i < n; i++) {
const int r = rev[i] >> shift;
if (i < r) std::swap(f[r], f[i]);
}
for (int s = 2, t = 1; s <= n; s <<= 1, t <<= 1) {
if (s == (1 << 18)) {
for (int i = 0; i < n; i++) f[i] %= mod;
}
static u64 val[1 << LIM_N - 1];
const u64 vI = rG[ctz(s)];
for (int i = val[0] = 1; i < t; i++) {
val[i] = val[i - 1] * vI % mod;
}
for (int x = 0; x < n; x += s) {
for (int y = 0; y < t; y++) {
const u64 vl = f[x | y], vr = f[x | y | t] * val[y] % mod;
f[x | y] = vl + vr, f[x | y | t] = vl - vr + mod;
}
}
}
for (int i = 0; i < n; i++) f[i] %= mod;
return f;
}
inline poly idft(poly f, int n) {
const int shift = LIM_N - ctz(n);
for (int i = 0; i < n; i++) {
const int r = rev[i] >> shift;
if (i < r) std::swap(f[r], f[i]);
}
for (int s = 2, t = 1; s <= n; s <<= 1, t <<= 1) {
if (s == (1 << 18)) {
for (int i = 0; i < n; i++) f[i] %= mod;
}
static u64 val[1 << LIM_N - 1];
const u64 vI = riG[ctz(s)];
for (int i = val[0] = 1; i < t; i++) {
val[i] = val[i - 1] * vI % mod;
}
for (int x = 0; x < n; x += s) {
for (int y = 0; y < t; y++) {
const u64 vl = f[x | y], vr = f[x | y | t] * val[y] % mod;
f[x | y] = vl + vr, f[x | y | t] = vl - vr + mod;
}
}
}
const u64 invn = qpow(n, mod - 2);
for (int i = 0; i < n; i++) f[i] = f[i] * invn % mod;
return f;
}
#undef ctz
inline poly cut(poly f, int n) { f.resize(n, 0); return f; }
inline poly Der(const poly &f) {
poly g(f.size() - 1);
for (int i = 0; i < g.size(); i++) {
g[i] = f[i + 1] * (i + 1) % mod;
}
return g;
}
inline poly Int(const poly &f) {
poly g(f.size() + 1); g[0] = 0;
for (int i = 0; i < f.size(); i++) {
g[i + 1] = f[i] * invi[i + 1] % mod;
}
return g;
}
inline poly operator + (const poly &f, const poly &g) {
poly h = f;
if (g.size() > f.size()) h = cut(h, g.size());
for (int i = 0; i < g.size(); i++) add(h[i], g[i]);
return h;
}
inline poly operator - (const poly &f, const poly &g) {
poly h = f;
if (g.size() > f.size()) h = cut(h, g.size());
for (int i = 0; i < g.size(); i++) add(h[i], mod - g[i]);
return h;
}
inline poly operator * (poly f, const u64 &v) {
for (u64 &x : f) x = x * v % mod;
return f;
}
inline poly operator * (poly f, poly g) {
const int n = f.size() + g.size() - 1;
int m = 1;
while (m < n) m <<= 1;
f = cut(f, m), g = cut(g, m);
f = dft(f, m), g = dft(g, m);
for (int i = 0; i < m; i++) f[i] = f[i] * g[i] % mod;
f = idft(f, m);
return f;
}
inline poly inv(const poly &h) {
assert(h[0] != 0);
const int n = h.size();
poly f{qpow(h[0], mod - 2)};
for (int m = 1; m < n; ) {
m *= 2;
f = cut(f * (poly{2} - cut(f * cut(h, m), m)), m);
}
return cut(f, n);
}
inline poly sqrt(const poly &h) {
assert(h[0] == 1);
const int n = h.size();
poly f{1};
for (int m = 1; m < n; ) {
m *= 2;
f = cut((f * f + cut(h, m)) * inv(cut(f, m)), m) * inv2;
}
f = cut(f, n);
return f;
}
inline std::pair<poly, poly> div(const poly &f, const poly &g) {
const int n = f.size(), m = g.size();
poly rf(f.rbegin(), f.rend()), rg(g.rbegin(), g.rend());
poly q = cut(rf * inv(cut(rg, n)), n - m + 1);
std::reverse(q.begin(), q.end());
poly r = cut(f - g * q, m - 1);
return std::make_pair(q, r);
}
inline poly ln(const poly &f) {
assert(f[0] == 1);
const int n = f.size();
return cut(Int(Der(f) * inv(f)), n);
}
inline poly exp(const poly &h) {
assert(h[0] == 0);
const int n = h.size();
poly f{1};
for (int m = 1; m < n; ) {
m *= 2;
f = cut(f * (poly{1} - ln(cut(f, m)) + cut(h, m)), m);
}
return cut(f, n);
}
inline poly pow_(const poly &f, const u64 &k1, const u64 &k2) {
assert(f[0] != 0);
const u64 f0 = f[0], invf0 = qpow(f0, mod - 2);
return exp(ln(f * invf0) * k1) * qpow(f0, k2);
}
inline poly pow(const poly &f, const u64 &k1, const u64 &k2, const u64 &k3) {
const int n = f.size();
if (k3 == 0) {
poly g(n, 0); g[0] = 1;
return g;
}
if (std::count(f.begin(), f.end(), 0) == n) {
return poly(n, 0);
}
if (f[0] == 0 && k3 > n) {
return poly(n, 0);
}
int index = -1;
for (int i = 0; i < n; i++) {
if (f[i]) {
index = i; break;
}
}
assert(index != -1);
poly g = pow_(poly(f.begin() + index, f.end()), k1, k2);
poly h(n, 0);
for (u64 i = k3 * index; i < n; i++) {
h[i] = g[i - k3 * index];
}
return h;
}
inline poly pow(const poly &f, const u64 &k) {
return pow(f, k % mod, k % (mod - 1), std::min(k, (u64)UINT32_MAX));
}
}
int main() {
ntt::init();
return 0;
}
后记
参考资料：
- 傅里叶变换(FFT)学习笔记 - command_block
- NTT与多项式全家桶 - command_block
- 算法学习笔记(13)：多项式牛顿迭代（求逆,开方,exp） - GhostLX
- 题解 P5245 多项式快速幂 - Constant
想放点私货失败了。
浙公网安备 33010602011771号