多项式入门到进阶超详解
前言：本文按照笔者自身理解撰写，力求通透易懂，适合萌新入门和大佬复习，如有不严谨处请多担待，有问题处欢迎反馈。最后如果您感觉有收获，还请您点个赞，谢谢啦。参考资料（同时也是笔者学会这些算法的出处）
https://oi-wiki.org/math/poly/fft/
https://oi-wiki.org/math/poly/ntt/
https://oi-wiki.org/math/poly/elementary-func/
复数：
我们都知道 \(x ^ 2 = 1\) 的解是：\(x_1 = 1,x_2 = -1\) 那么方程 \(x ^ 4 = 1\) 的解是什么呢？
首先一个 \(N\) 次方程如果有 \(N\) 个解，所以这里应该有四个解。
因为 \(x ^ 4 = (x ^ 2) ^ 2\) 所以 \(x ^ 2 = 1\) 或 \(x ^ 2 = -1\) 这时我们会发现 \(x = \sqrt{-1}\) 我们称这样的 \(x\) 叫做 \(i\)，同理 \(i\) 也满足 \(i ^ 2 = -1\)。
这里引入复数：对于任意的数，都可以表示成 \(a + bi\) 这里\(a,b\in R\) 并且实数是完全被复数包含在内的。
我们知道，对于任意实数 \(x\in R\) 其在数轴上都能被唯一地表示出来。那么对于复数，我们应该用什么样的坐标系，才能让他们都能被唯一地表示出来呢？
考虑每个 \(a + bi\) 我们将其常数项和 \(i\) 的系数组成一个二元组，发现对于每种二元组都能唯一地对应一个复数，或者说，每个这样的二元组都和其对应的复数组成双射关系，那么我们就可以把它们放到二维坐标系中。
一些概念：
[*]我们把复数 \(a + bi\) 对应的坐标 \((a,b)\) 所在的坐标系中的横轴表示实轴，纵轴表示虚轴。
[*]如果我们把原点 \(O\) 和一个复数 \(a + bi\) 对应的坐标 \((a,b)\) 连边，那么这条边的长度就称作这个复数的模长。
[*]如果我们把原点 \(O\) 和一个复数 \(a + bi\) 对应的坐标 \((a,b)\) 连边，那么这条边和实轴的正方向所夹的角称作幅角。
由于每个复数所对应的二元组在二维坐标系中，我们可以将它和向量的知识联系在一起。
复数运算：
[*]加法运算，设我们有两个复数 \(a + bi\) 和 \(c + di\) 那么运算法则就是对应位置（即：常数相加，\(i\) 的系数相加）相加，也就是 \(a + bi + c + di \rightarrowa + c + bi + di \rightarrow a + c + (b + d)i\)
[*]减法运算，其实就是加法运算的逆，可以把 \(a + bi - (c + di)\rightarrow a + bi + (-c - di)\) 即 \(a + bi - (c + di) \rightarrow a - c + (b - d)i\)
[*]乘法运算，其实就是 \((a + bi)\times (c + di) \rightarrow ac + bdi^2 + bci + adi\) 因为 \(i ^ 2 = -1\)，所以原式就是：\((a + bi)\times (c + di) \rightarrow ac - bd + (bc + ad)i\)
[*]除法运算，其实就是 \(\frac{a + bi}{c + di}\) 可以分子分母同乘 \(c - di\) 从而分子变成平方差，即 \(\frac{a + bi}{c + di} \times \frac{c - di}{c - di} \rightarrow \frac{ac - adi + bci -bdi^2}{c ^ 2 - d^2i^2}\rightarrow \frac{ac + bd + (bc - ad) i}{c ^ 2 + d ^ 2}\)
[*]复数的共轭：实部相同，虚部相反(根据实轴对称)。
[*]复数相乘的时候模长相乘，幅角相加。
单位圆：
和三角函数一样，就是以 \(O\) 为圆心，半径为一个单位长度的圆。
不难发现，单位圆上的每一个点模长都是 \(1\)
为了方便书写，我们下文称 \(|x|\) 为 \(x\) 的模长。
快速傅里叶变换思想：
引理： 一个 \(N\) 次多项式可以由 \(N + 1\) 个 \(x\) 互不相同的点值唯一确定出来，比如一次函数需要两个点，二次函数需要三个点。
一些符号的说明： 设当前有两个多项式 \(f,g\)，我们称 \(fg\) 表示多项式 \(f\) 和 \(g\) 相乘，即，若 \(h = fg\) 则 \(h_i = \sum_{j + k = i} f_j\times g_k\)，定义 \(deg_f\) 表示多项式 \(f\) 次数。
首先暴力做法就是对于 \(0\leq i\leq deg_f,0\leq j\leq deg_g, h_{i + j}\rightarrow h_{i + j} + f_i\times g_j\)
时间复杂度是 \(O(n^2)\) 的。
我们设多项式 \(f(x) = a_0 + a_1x + a_2x^2 + ... a_nx^n\) 中系数为 \(\) 我们发现直接对系数进行合并是难以优化的，但是如果转成点值呢？假设我们已经把 \(f,g\) 分别转成点值序列 \(A\) 和 \(B\) 那么点值相乘是 \(O(n)\) 的，即 \(C = A\times B,\forall i\in ,C_i = A_i\times B_i\) 最后我们只需要再进行一次多项式转点值的逆，也就是点值转多项式就好了，所以我们把重点放到了如何知道一个多项式 \(f\) 求出 \(deg_f + 1\) 个点的点值序列。
离散傅里叶变换 DFT
既然我们打算要把一个多项式 \(f\) 转成一个点值序列 \(A\) 那么我们就需要任意找 \(deg_f + deg_g + 1\) 个点分别带入两个多项式，因为 \(|fg| = |f| + |g| - 1\) 而 \(|f| = deg_f + 1\) 其中 \(|f|\) 表示多项式 \(f\) 有多少个幂， \(0\) 次幂也算幂。
但是想法是好的，显示是很残酷的，因为多项式多点求值需要使用 \(FFT\) 实现，但是我们研究的是 \(FFT\) 肯定是不能这样的，因此我们应该找一些特殊点。
伟大的傅里叶告诉我们，要从单位根上来考虑这个问题。
具体的，我们把单位圆均分成 \(N\) 份。
我们发现，这 \(N\) 个根也是 \(x ^ N = 1\) 的解，后文称作单位根，我们设实轴上的边为 \(\omega_{n}^0\) 然后逆时针的过程中第一个遇到的红线就是 \(\omega_{n}^1\) 同理，第 \(i\) 个遇到的红线就是 \(\omega_{n}^i\) 容易发现，由于我们有 \(n\) 个单位根，所以也就是有 \(\omega_{n}^0 .... \omega_{n}^{n - 1}\) 下面我们重点来讨论一下这些单位根的性质。
单位根性质：
[*]\(\omega_n^{0,1.....n - 1}\) 并不是几个独立的点，而是代表一个集合，因为如果在某一个点，它逆时针再走 \(2\pi\) 步仍然可以到这个点，由此也能推出单位圆中的周期性 \(T = 2\pi\) 或者说 \(w_n^i = w_n^{i\bmod n},i\in Z\)
[*]\(\omega_n^j \times \omega_n^i = \omega_n^{i + j}\) 这个就是刚刚复数运算中的复数相乘，幅角相加。
[*]\((\omega_n^i)^j = \omega_n^{ij}\) 这个其实就是 \(j\) 个 \(w_n^i\) 相乘和 \(2\) 一样的分析。
[*]\(\omega_{2n}^{2i} = \omega_n^i\) 类比分蛋糕的时候，就是原来是分了 \(2n\) 块然后选了 \(2i\) 块，现在是分了 \(n\) 块，选了 \(i\) 块，不难发现是一样的，注意：这里的“分”指的是均分。
[*]如果 \(n\) 是偶数，那么 \(\omega_n^i = -\omega_n^{i + \frac{n}{2}}\) 因为 \(n\) 是偶数，所以 \(w_n^{\frac{n}{2}}\) 就相当于逆时针走了 \(180\) 度，自然就是变成了原来中心对称的位置了。
准备好了吗，即将开始最重要的部分！
考虑当前有个 \(n\) 次多项式 \(f(x) = a_0 + a_1x + a_2x^2 + ... a_nx^n\) 并且这里 \(n - 1\) 是偶数。
我们把他按照奇偶性分开：
令 \(g(x) = a_0 + a_2x^2 + a_4x^4 + ... a_{n - 2}x^{n - 2}\)
令 \(h(x) = a_1x + a_3x^3 + a_5x^5 + ... a_{n - 1}x^{n - 1}\)
显然 \(f(x) = g(x) + h(x)\) 这里 + 是按位相加。
再变：
\(g(x)\rightarrow a_0 + a_2x + a_4x^2 + a_{n - 2}x^{\frac{n}{2} - 1}\)
\(h(x)\rightarrow a_1 + a_3x + a_5x^2 + a_{n - 1}x^{\frac{n}{2} - 1}\)
所以现在 \(f(x) = g(x ^ 2) + xh(x ^ 2)\)，我们就可以带值找找规律了。
设当前带入 \(\omega_n^k,k < \frac{n}{2}\) 分两种情况带入：
[*]带入 \(\omega_{n}^k\) 此时 \(f(\omega_{n}^k) = g((\omega_{n}^k)^2) + \omega_{n}^kh((\omega_{n}^k)^2)\)
[*]\(f(\omega_{n}^k) = g((\omega_{n}^k)^2) + \omega_{n}^kh((\omega_{n}^k)^2)\rightarrow f(\omega_{n}^k) = g(\omega_{n}^{2k}) + \omega_{n}^kh(\omega_{n}^{2k})\rightarrow f(\omega_{n}^k) = g(\omega_{\frac{n}{2}}^k) + \omega_{n}^kh(\omega_{\frac{n}{2}}^k)\)。
[*]由此得出来了，对于 \(k < \frac{n}{2}\)，\(f(\omega_n^k)\) 的点值就是 \(g(\omega_{\frac{n}{2}}^k) + \omega_n^kh(\omega_{\frac{n}{2}} ^ k)\) 也就是说，而 \(g,h\) 都是子问题，且范围缩小了一倍，所以我们每次把原问题的范围分成两半来分别解决然后合并，这是典型的分治，总的层数是 \(O(\log n)\) 级别的。
[*]带入 \(\omega_{n}^{k + \frac{n}{2}}\) 此时 \(f(\omega_{n}^{k + \frac{n}{2}}) = g((\omega_{n}^{k + \frac{n}{2}}) ^ 2) + \omega_n^{k + \frac{n}{2}}h((\omega_n^{k + \frac{n}{2}}) ^ 2)\)
[*]因为 \((\omega_{n}^j)^k = \omega_n^{jk}\)，原式得：\(f(\omega_{n}^{k + \frac{n}{2}}) = g(\omega_{n}^{2k + n}) + \omega_n^{k + \frac{n}{2}}h(\omega_n^{2k + n})\)
[*]因为 \(\omega_n^j = \omega_n^{j\bmod n}\)，原式得：\(f(\omega_{n}^{k + \frac{n}{2}}) = g(\omega_{n}^{2k}) + \omega_n^{k + \frac{n}{2}}h(\omega_n^{2k})\)
[*]因为 \(\omega_{2n}^{2j} = \omega_n^j\)，原式得：\(f(\omega_n^{k + \frac{n}{2}}) = g(\omega_{\frac{n}{2}} ^ k) + \omega_n^{k + \frac{n}{2}} h(\omega_{\frac{n}{2}}^k)\)
[*]这里为了方便，可以根据 \(\omega_n^{k + \frac{n}{2}} = -\omega_n^k\)，原式得：\(f(\omega_n^{k + \frac{n}{2}}) = g(\omega_{\frac{n}{2}} ^ k) - \omega_n^{k} h(\omega_{\frac{n}{2}}^k)\)
因此我们就得到了关于计算 \(f(\omega_n^i)\) 的式子，也就是每次把原序列折半分成奇数偶数两半，然后分别计算子问题的 \(f(\omega_n^i)\) 最后合并，不难发现我们转换后的式子恰好满足 \(n\rightarrow \frac{n}{2}\) 且 \(k < \frac{n}{2}\) 所以成功划分成了子问题！
这里有些细节：
如果递归到某一层， \(n = 1\) 那么直接 return，这个是因为递归到 \(n = 1\) 的时候是计算 \(f(w_1^0)\) 由于 \(w_1 ^ 0\) 其实就是 \(1 + 0i = 1\) 所以 $f(w_1 ^ 0) = $ 分治到 \(i\) 时 \(i\) 的系数，所以可以直接分治的时候记录系数，然后把系数当做 \(f(w_1 ^ 0)\) 来计算。
时间复杂度不难发现是 \(T(n) = 2T(\frac{n}{2}) + O(n) = O(n\log n)\)
离散逆傅里叶变换 IDFT
假设我们已经进行了一次正变换 DFT，得到的 \(N\) 个点值是序列 \(\{y\},\forall i\in\) 对于 \(y_i = f(\omega_n^i)\) 设原来的函数是 \(f(x) = a0 + a1x + a2x ^ 2 + a3x ^ 3 + ... a_{n - 1} x ^ {n - 1}\) 我们现在来构造一个新的多项式 \(A\) 满足：
\(A(x) = \sum_{0\leq i< n} y_ix^i\)
我们再设 \(b\) 序列，其中 \(b_i = \omega_n^{-i}\) 那么我们来看看把 \(\{b\}\) 序列带入 \(A(x)\) 看看会产生什么样的点值序列。
\(A(b_k) = \sum_{i = 0}^{n - 1} y_i{b_k}^i\)
因为 \(y_i = f(\omega_n ^ i) =\sum_{j = 0}^{n - 1} a_j({\omega_n ^ i}) ^ j\)
所以原式得：\(A(b_k) = \sum_{i = 0}^{n - 1} {b_k} ^ i \sum_{j = 0}^{ n- 1} a_j(\omega_n^{i})^j\)
因为： \((\omega_n^i)^j = \omega_n ^ {ij}\)
所以原式得：\(A(b_k) = \sum_{i = 0}^{n - 1} {b_k} ^ i \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}\)
因为 \(b_k = \omega_n^{-k}\)
所以原式得：\(A(b_k) = \sum_{i = 0}^{n - 1} {\omega_n^{-k}} ^ i \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}\rightarrow \sum_{i = 0}^{n - 1} {\omega_n^{-ik}} \sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}\)
交换求和：
\(A(b_k) = \sum_{i = 0}^{n - 1}\sum_{j = 0}^{ n- 1} {\omega_n^{-ik}}a_j\omega_n^{ij}\)
因为：\(\omega_n^i\omega_n^j = \omega_n^{i + j}\)
所以原式得：\(A(b_k) = \sum_{i = 0}^{n - 1}\sum_{j = 0}^{ n- 1} a_j\omega_n^{ij}\omega_n^{-ik}\rightarrow \sum_{i = 0}^{n - 1}\sum_{j = 0}^{ n- 1} a_j\omega_n^{i(j - k)}\)
因为 \((\omega_n^i)^j = \omega_n^{ij}\)
所以原式得：\(A(b_k) = \sum_{i = 0}^{n - 1}\sum_{j = 0}^{ n- 1} a_j(\omega_n^{j -k}) ^ i\)
交换求和得：
\(A(b_k) = \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1} (\omega_n^{j - k})^i\)
此时记 \(S(\omega_n^a) = \sum_{i = 0}^{n - 1} (\omega_n^{a}) ^ i\)
[*]如果 \(a \bmod n = 0\) 那么显然 \(\omega_n^a = 1\) 所以 \(\sum_{i = 0}^{n - 1}(\omega_n^a) ^ i = n \times 1 = n\)
[*]否则考虑相邻两项作差：
[*]等式两边同乘 \(\omega_n^a\)
[*]原式得:\(\omega_n ^ a S(\omega_n^a)=\sum_{i = 1}^{n}(\omega_n^{a}) ^ i\)
[*]两者作差得：\(\omega_n^aS(\omega_n^a) - S(\omega_n ^ a) = \sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i\)
[*]左边提取公因式得：\((\omega_n^a - 1)S(\omega_n^a)= \sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i\)
[*]由于 \(a\bmod n\ne 0\) 所以 \(\omega_n^a\ne 1\) 所以 \((\omega_n^a - 1)\ne 0\) 所以可以等式两边同时除以 \((\omega_n^a - 1)\)
[*]因此，原式得：\(S(\omega_n^a) = \frac{\sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i}{\omega_n^a - 1}\)
[*]不难发现 \(\sum_{i = 1}^n (\omega_n^a)^i - \sum_{i = 0}^{n - 1}(\omega_n^a)^i = (\omega_n^a)^n - (\omega_n^a)^0\)
[*]前者可以写成 \((\omega_n^n)^a = (\omega_n^0)^n = 1\) 后者是 \(0\) 次方也是 \(1\)
[*]所以可以得到对于 \(a\bmod n \ne 0\) 的所有 \(a\) 满足 \(S(\omega_n^a) = 0\)
所以最后可以归纳成：对于 \(a < n\) 的所有 \(\omega_n^a\) 满足：也就是说
$ S\left(\omega_n^a\right)=
\begin{cases}
n,&a=0\
0,&a\neq 0
\end{cases}
$
好的，现在我们把这个公式带到 \(A(b_k) = \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1} (\omega_n^{j - k})^i\)
原式得：\(\sum_{j = 0}^{n - 1} a_j S(\omega_n^{j - k})\)
显然这个式子，只有在 \(j = k\) 的时候有值，且值为 \(n\) 所以这个式子的值就是 \(na_k\)
也就是说，给定点 \(b_i = \omega_n ^ {-i}\) 它的点值序列表示就是 \(\{(b_0,a_0n),(b_1,a_1n),....,(b_{n - 1},a_{n - 1}n)\}\)
所以我们就得到了如何在求出 \(f,g\) 的点值序列后如何还原出来新的多项式的系数，当然，因为最后还原出来的点值是 \(a_in\) 所以要除以 \(n\)。
不难发现，我们这个过程和 DFT 几乎一模一样，无非就是将 \(n\) 个单位根 \(\omega_n^i\rightarrow \omega_n^{-i}\) 这两个式子是关于实轴（x轴）对称的，所以只需要把纵坐标乘上一个 \(-1\) 即可。
但是我们发现直接写分治，我们空间复杂度应该也是 \(O(n \log n)\) 的，因为需要对每个点开一个 \(O(len)\) 级别的数组，所以要引进一种空间上的优化：蝶形优化。
蝶形优化
对于当前分治区间 \(\) 设当前长度是 \(n_i = r_i - l_i + 1\)
我们要计算 \(f(\omega_{n_i}^{0,1,...n_i - 1})\) 需要 \(g(\omega_{\frac{n_i}{2}} ^ {0,1,2,...\frac{n_i}{2} - 1})\) 和 \(h(\omega_{\frac{n_i}{2}} ^ {0,1,2,...\frac{n_i}{2} - 1})\)
我们可以在递归开始前，就把 \(g,h\) 给放到 \(f\) 的第 \(\) 的位置和 \([\frac{n_i}{2},n_i - 1]\) 的位置上，这个位置是个指针，此时修改 \(g\) 和 \(h\) 也会在对应的 \(f\) 上进行修改，（类似长剖一样）所以在我们更新 \(f\) 的时候不能直接拿 \(g\) 和 \(h\) 来改 \(f\) 而是新开一个辅助数组 \(G\) 最后再把 \(G\) 赋值给 \(f\) 不然直接改 \(f\) 的值，对应的 \(g\) 和 \(h\) 的值也会进行修改，就不对了。这样空间就被优化成 \(O(n)\) 级别了。
最后，一些实现上的小细节：
[*]我们发现当 \(n = 1\) 的时候 \(\omega_n^{n - 1} = 1\) 的，所以可以直接把一开始的点值放进去，因为此时点值的 \(f\) 值就是自己不会变。所以一开始的时候，就直接输入 \(A,B\) 的实部就可以了。
[*]一定要注意，\(f(i) = a + bi\) 中 \(a + bi\) 中的 \((a,b)\) 是个整体，代表一个复数，而不是当 \(x = a\) 的时候对应的值是 \(b\)。
[*]分治的时候，一定不要直接修改 \(f\) 的值，因为使用蝶形优化后，\(g\) 和 \(h\) 和 \(f\) 公用一个数组。
[*]因为我们刚才说的，每次分治的时候都需要能均分奇偶数，所以要求所有大于 \(1\) 的分治区间长度是 \(2\) 的倍数，自然想到构造一个 \(2 ^ k\) 的数列，如果一开始不够就补 \(0\) 即可，然后找到最大的 \(k\) 满足 \(2^k \ge n + m + 1\)
[*]由于我们是进行的实数运算，会有精度问题，所以最后要输出四舍五入后的值。
代码：
#include #include #include constexpr double Pi = acos(-1);constexpr int N = 4e6 + 5;#define int long longstruct complex{ double x,y; inline complex operator + (const complex & a){ return {a.x + x,a.y + y}; } inline complex operator - (const complex & a){ return {x - a.x,y - a.y}; } inline complex operator * (const complex & a){ // (x + yi) * (a.x + a.yi) // a.x * x + xa.yi + a.xyi + a.yi*yi // a.x * x + xa.yi + a.xyi - a.y*y return {a.x * x - a.y * y,x * a.y + a.x * y}; }}A,B,g;int n , m;inline void fft(complex * f,int len,int rev){ if(len == 1) return; for(int i = 0; i < len; ++ i){ g = f; } complex * ls = f,* rs = f + len / 2; for(int i = 0; i < len / 2; ++ i){ ls = g.x; } for(int i = 0; i > B.x; } int len = 1; for(int k = n + m; lenm; for(int i = 0; i > a.x ; // 读入实部 } for(int i = 0; i > b.x; } for(int k = n + m; len >= 1,(x *= x) %= mod){ if(y & 1) (ans *= x) %= mod; } return ans;}inline void NTT(ll f[],int len,int rev){ reverse(f,len); for(int L = 2; L > m; for(int i = 0; i > a; } for(int i = 0; i > b; } for(int k = n + m; len1; if(i & 1){ R |= len >> 1; } } for(int i = 0; i < len; ++ i){ if(i < R) std::swap(f,f]); }}inline void NTT(ll *f,int len,int rev){ reverse(f,len); for(int L = 2; L >= 1,(x *= x) %= mod){ if(y & 1) (ans *= x) %= mod; } return ans;}struct poly{ int R,c,d; inline void reverse(int f[],int n){ R = 0; for(int i = 1; i < n; ++ i){ R = R >> 1; if(i & 1){ R |= n >> 1; } } for(int i = 0; i < n; ++ i){ if(i < R){ std::swap(f,f]); } } } inline void NTT(int f[],int n,int rev){ reverse(f,n); for(int L = 2; L
页:
[1]