多项式乘法
多项式算法作为 OI 中较为困难的部分，其包含的知识是比较丰富的。本文将介绍多项式算法的基础——多项式乘法。
1 快速傅里叶变换 FFT
1.1 引入
首先我们知道一个多项式 \(f(x)\) 可以写作：
那么我们如果想对两个多项式做乘法，会有如下结果：
也就是说，得到的新多项式的系数可以写成：
这种形式被称为卷积，其特征是相乘的两个数下标为定值。不难发现，如果直接暴力计算，对于两个 \(n\) 次多项式来讲，其复杂度是 \(O(n^2)\) 的，并不优秀。
为了实现以后的更多多项式操作，我们需要降低多项式乘法的复杂度，于是就有了快速傅里叶变换（Fast Fourier Transform，FFT）来帮助我们求解这个问题。
1.2 前置知识
1.2.1 系数表示与点值表示
FFT 的基础是系数表示和点值表示。系数表示不必多说，对于一个多项式 \(f(x)=\sum\limits_{i=0}^n a_ix^i\)，我们将 \(a\) 数组可以看作一个 \(n+1\) 维向量 \(\overrightarrow{a}=(a_0,a_1,\cdots,a_n)\)，那么该多项式的系数表示就是 \(\overrightarrow{a}\)。
重要的部分是点值表示。在拉格朗日插值中我们提到过这样一个定理：对于一个 \(n\) 次多项式，只需要给定 \(n+1\) 个点，就可以唯一确定这个多项式。因此我们取 \(n+1\) 个不同的 \(x_0,x_1,\cdots x_n\) 代入 \(f(x)\) 可以得到 \((x_0,f(x_0)),(x_1,f(x_1)),\cdots,(x_n,f(x_n))\) 这些点，它们就可以唯一代表这个多项式。
这有什么用呢？考虑两个多项式 \(f(x)\) 与 \(g(x)\)，它们卷积的系数表示上文已经说过，但是从点值表示来看，它们的卷积结果就是 \((x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),\cdots,(x_n,f(x_n)g(x_n))\)。不难发现此时我们只需要 \(O(n)\) 的复杂度就可以直接求出两个多项式卷积的结果了。现在的问题就是如何在多项式的系数表示和点值表示之间进行快速的转换。
1.2.2 复数
复数是形如 \(a+bi\) 的数，其中 \(i=\sqrt{-1}\)，即复数单位。其中 \(a\) 称为实部，\(b\) 称为虚部。
然后可以定义复数的加、减、乘运算如下：
- \((a+bi)+(c+di)=(a+c)+(b+d)i\)。
- \((a+bi)-(c+di)=(a-c)+(b-d)i\)。
- \((a+bi)(c+di)=ac+adi+bci-bd=(ac-bd)+(ad+bc)i\)。
接下来我们引入复平面，以实部作为 \(x\) 轴，虚部作为 \(y\) 轴建立平面直角坐标系，于是任何一个复数都可以表示为复平面上的一个点。
1.2.3 单位根
对于一个 \(n\) 次方程 \(x^n=1\) 来讲，它应该会有 \(n\) 个根。我们将这个方程的所有根对应到复平面上，会发现它们呈现如下形式：以原点为圆心，做一个半径为 \(1\) 的圆，称为单位圆；则这些根会 \(n\) 等分单位圆，且始终有一个根为 \(1\)。
我们称这 \(n\) 个根为 \(n\) 次单位根，同时将从 \(1\) 开始逆时针找到的第一个根记作 \(\omega_n\)。于是不难发现，剩下的所有单位根都可以用 \(\omega_n\) 的幂次表示，即 \(n\) 次单位根分别为 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\)。
根据定义，可以利用三角函数简单求出 \(\omega_n\) 的表达式：
然后我们还需要知道单位根的若干重要性质，这在 FFT 中有很大用处：
- \(\omega_n^k=\omega_n^{k+n}\)。
- \(\omega_n^k=\omega_{2n}^{2k}\)。
- \(\omega_{2n}^{n+k}=-\omega_{2n}^k\)。
证明都不困难，结合复平面上的几何意义以及适量代数推导都可以简单得出。
1.3 离散傅立叶变换 DFT
离散傅里叶变换 DFT 实际上描述的就是将多项式从系数表示转化成点值表示的一种方法。例如对于多项式 \(f(x)=\sum\limits_{i=0}^{n-1}a_ix^i\)，我们将 \(n\) 个不同的 \(n\) 次单位根代入多项式 $f(x) $ 中，就可以得到其对应的点值表达。不过这样暴力做的复杂度仍然是 \(O(n^2)\) 的。而 FFT 就提供了一种快速进行 DFT 的算法。
考虑一个 \(2n-1\) 次多项式 \(f(x)\)，如果我们分别设：
那么我们显然有 \(f(x)=g(x^2)+xh(x^2)\)，此时代入 \(x=\omega_{2n}^k\)，可得：
同理，将 \(x=\omega_{2n}^{n+k}\) 代入，可得：
于是当我们求出 \(g(\omega_n^k)\) 和 \(h(\omega_n^k)\) 之后，可以同时求出 \(f(\omega_{2n}^k)\) 和 \(f(\omega_{2n}^{k+n})\)。此时我们将 \(f\) 拆成了两个 \(n-1\) 次多项式，于是可以直接递归对 \(g,h\) 再做 DFT 即可。复杂度显然是 \(O(n\log n)\)。
实现时为了方便，我们需要将多项式的项数补到一个 \(2\) 的整数次幂，方便我们分治处理。
当然实际运用中一般不会写递归分治写法，原因在于其常数过大，我们一般采用倍增法实现。关于倍增实现有两个基本要点：位逆序置换和蝶形运算优化。
考虑我们拆分多项式的过程，设有一个 \(8\) 项多项式，其系数编号会如下变化：
- 初始下标为 \(\{0,1,2,3,4,5,6,7\}\)。
- 第一次分治后下标为 \(\{0,2,4,6\}\{1,3,5,7\}\)。
- 第二次分治后下标为 \(\{0,4\}\{2,6\}\{1,5\}\{3,7\}\)。
- 第三次分治后下标为 \(\{0\}\{4\}\{2\}\{6\}\{1\}\{5\}\{3\}\{7\}\)。
不妨写出最开始和最后系数下标的二进制表示，会发现如下情况：
- 初始下标为 \(000,001,010,011,100,101,110,111\)。
- 结束下标为 \({000,100,010,110,001,101,011,111}\)。
对比后不难发现，每一个下标都只是翻转了其二进制表示。这就是位逆序置换。我们可以提前 \(O(n)\) 预处理出每一个数字位逆序置换后的结果， 然后将每一对要交换的数交换即可。
接下来是蝶形运算优化，我们会发现在拆分系数之后，\(g(\omega_n^k)\) 存储在第 \(k\) 项，而 \(h({\omega_{n}^{k}})\) 存储在第 \(k+n\) 项。在合并之后，\(f(\omega_{2n}^k)\) 存储在第 \(k\) 项，\(f(\omega_{2n}^{n+k})\) 存储在第 \(k+n\) 项。于是我们只需要对第 \(k\) 项和第 \(k+n\) 项进行覆写即可。
1.4 逆离散傅立叶变换 IDFT
和上面的 DFT 类似，IDFT 就是将点值表示转化为系数表示的变换。不过我们不必实现一个新的函数，因为实际上 IDFT 可以转化为 DFT 来实现。
上述 DFT 的过程可以看成是对系数表示的向量左乘了一个矩阵，得到了点值表示的向量，即下式：
将中间的大矩阵记作 \(\mathbf{V}\)，现在构造一个新矩阵 \(\mathbf{D}\)，如下：
即 \(d_{i,j}=(\omega_n^{-i})^j\)。然后让我们令 \(\mathbf{E}=\mathbf{D}\cdot\mathbf{V}\)，可以得到：
接下来不难发现：
- 当 \(i=j\) 时，显然 \(e_{i,j}=n\)。
- 当 \(i\ne j\) 时，运用等比数列求和公式可知 \(e_{i,j}=\dfrac{(\omega_n^{j-i})^n-1}{\omega_n^{j-i}-1}=\dfrac{(\omega_n^n)^{j-i}-1}{\omega_n^{j-i}-1}=0\)。
所以说 \(\mathbf{E}\) 矩阵只有主对角线上的数字为 \(n\)，其它全部为 \(0\)。也就是说，\(\mathbf{E}\) 矩阵是单位矩阵的 \(n\) 倍。因此 \(\mathbf{D}\) 的 \(\dfrac 1n\) 就是 \(\mathbf{V}\) 的逆矩阵。于是在等式两边同时左乘上 \(\dfrac 1n\mathbf{D}\) 即可得到点值表示向量化为系数表示向量的方法了。
由于上面我们实际上 \(\mathbf{D}\) 只是将原先的 \(\omega_n^k\) 换成了 \(\omega_n^{-k}\)，所以只需要将 DFT 函数中的单位根换成 \(\omega_n^{-k}\)（实际上只需要对虚部取相反数即可，证明并不困难），最后将结果除以 \(n\) 即可实现 IDFT。
1.5 常数优化及代码
在上述过程中，我们对两个多项式分别作了一次 DFT，乘起来之后又作了一次 IDFT，总共三次 FFT 操作。实际上我们可以通过一些操作使得只需要进行两次 FFT 操作。
具体的，我们将第二个多项式的系数放到第一个多项式的虚部上，然后让第一个多项式平方。根据 \((a+bi)^2=(a^2+b^2)+2abi\)，平方后我们只需要将虚部除以 \(2\) 即可得到正确的答案。
不过这种方法也有弊端，如果两个相乘的多项式值域相差过大，那么由于有平方项的存在这种做法会严重掉精度，然后导致 WA。解决方法就是将两个多项式的系数先乘到一个差不多一致的值域范围内，求答案的时候再除回去即可。不过这样的话除法还有巨大的常数，所以使用的时候需要考虑一下。
模板题：【模板】多项式乘法（FFT），采用朴素三次 FFT 实现代码如下：
#include <bits/stdc++.h>
using namespace std;
const int Maxn = (1 << 21) + 1;
const int Inf = 2e9;
const double Pi = acos(-1);
struct comp {//手写复数类
double a, b;
comp operator + (comp y) {return {a + y.a, b + y.b};}
comp operator - (comp y) {return {a - y.a, b - y.b};}
comp operator * (comp y) {return {a * y.a - b * y.b, a * y.b + b * y.a};}
comp operator / (int y) {return {a / y, b / y};}
};
int r[Maxn];
struct Poly {
int n;//n 次多项式
vector <comp> a;
void reset(int len) {n = len, a.resize(len + 1);}//调整多项式长度
comp& operator [](int x) {return a[x];}
void FFT(int len, int typ) {
//处理 len 项多项式
reset(len - 1);
for(int i = 0; i < len; i++) {
if(i < r[i]) swap(a[i], a[r[i]]);
}
for(int h = 1; h < len; h <<= 1) {
//枚举小区间长度 h，当前处理长度为 2h 区间
comp cur = {cos(Pi / h), typ * sin(Pi / h)};//单位根 w(2h)
for(int i = 0; i < len; i += (h << 1)) {//跳大区间
comp w = {1, 0};
for(int j = 0; j < h; j++, w = w * cur) {//遍历小区间
comp x = a[i + j], y = w * a[i + j + h];
a[i + j] = x + y;
a[i + j + h] = x - y;//蝶形运算
}
}
}
if(typ == -1) {//IDFT 的时候要除以项数 len
for(int i = 0; i < len; i++) a[i] = a[i] / len;
}
}
Poly operator * (Poly y) {
Poly x = *this, z;
int n = x.n + y.n, len = 1;
while(len <= n) len <<= 1;//将项数变成 2 的整数次幂
for(int i = 0; i < len; i++) {//预处理位逆序
r[i] = (r[i >> 1] >> 1) | ((i & 1) * (len >> 1));
}
x.FFT(len, 1), y.FFT(len, 1);
z.reset(len);
for(int i = 0; i < len; i++) z[i] = x[i] * y[i];
z.FFT(len, -1);
z.reset(n);
return z;
}
}F, G;
int n, m;
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> m;
F.reset(n), G.reset(m);
for(int i = 0; i <= n; i++) cin >> F[i].a;
for(int i = 0; i <= m; i++) cin >> G[i].a;
F = F * G;
for(int i = 0; i <= n + m; i++) cout << (int)(F[i].a + 0.5) << ' ';//精度问题
return 0;
}
1.6 例题
例 1 [BZOJ2194] 快速傅里叶之二
题意：求 \(C_k=\sum\limits_{i=k}^{n-1} A_i\times B_{i-k}\)。
容易发现此时后面的下标差值恒为定值 \(k\)，但是我们卷积要求是下标和为定值。
为了使下标和为定值，我们需要将里面的一个 \(i\) 改成 \(-i\)。于是不难想到如下方法：将数组 \(A\) 翻转，这样 \(A_i\) 就可以变为 \(A_{n-i-1}\) 了。此时做卷积可得 \(C'_{n-1-k}=\sum\limits_{i=k}^{n-1} A_{n-i-1}\times B_{i-k}=\sum\limits_{i=0}^{n-1-k} A_{n-i-k-1}\times B_{i}\)。此时直接对 \(A,B\) 做卷积，然后 \(C'_{n-1}\sim C'_{0}\) 就对应 \(C_0\sim C_{n-1}\) 了。
实际上这就告诉我们差卷积也是卷积的基本形式，也可以直接 FFT 求解
例 2 [Tyvj1953] Normal
首先难点实际上在于转化期望式子，转化不对你是做不出来的。我们考虑每一个点都会做一次分治中心，当点 \(i\) 做分治中心的时候，其贡献就是其子树内点的个数。根据期望线性性，不难得出答案就是下式：
考虑后面的概率怎么求，不难发现的是，如果 \(j\) 在 \(i\) 的分治子树内，那么 \(i\) 到 \(j\) 路径上的这些点中最早被删除的一定是 \(i\)。而这个概率就是 \(\tfrac1{dis(i,j)}\)。于是答案就是：
考虑怎样求这个，实际上就是求每种长度的路径有多少条，这个时候就要用点分治了。考虑怎样统计子树信息，由于子树合并的复杂度我不会证，所以考虑子树容斥。具体的，先求出当前分治中心内所有点到分治中心的距离，然后我们需要两两组合拼成一条新路径。设距离分治中心距离为 \(i\) 的点有 \(p_i\) 个，那么合并后长为 \(n\) 的路径数就应该有：
标准的不能再标准的卷积形式了，将 \(p\) 自己卷自己即可。然后进行子树容斥，把路径端点在同一子树内的方案数减掉，实际上就是对每个子树再像上面做一遍即可。复杂度 \(O(n\log^2 n)\)。
例 3 [BZOJ3771] Triple
题意： 给出 \(n\) 个数构成的序列，从中选出 \(1\) 个、\(2\) 个、\(3\) 个数，求所有可能的权值和及方案数。两种方案本质不同当且仅当选出的数的下标构成的集合不同。
其实这道题运用的思想和生成函数差不多。我们构造这样一个多项式 \(F\)：对于每一项 \(a_ix^i\)，\(a_i\) 表示序列中值为 \(i\) 的数的个数。这样的话我们将两个多项式做卷积，指数相加系数相乘，正好可以实现权值和、方案数相乘的效果。所以多项式乘法是可以用于解决背包问题的。
然后考虑怎么求，分类讨论即可。
- 对于选 \(1\) 个，直接累加到答案里即可。
- 对于选 \(2\) 个，先求出 \(F\times F\)，然后需要减去选两个重复下标的方案数，重新构造一个多项式 \(G\) 来表示这个方案数，答案即为 \(F\times F-G\)，注意最后要除以 \(2\)。
- 对于选 \(3\) 个，先求出 \(F\times F\times F\)，然后要减去选了两个重复下标以及三个重复下标的方案数。构造一个多项式 \(H\) 来表示选三个重复下标的方案数。对于前者先求出 \(F\times G\)，但是这样又多算了选三个重复下标的方案，所以再减去 \(H\)；并且由于这种选法会重复算三次，所以应该减去 \(3\times (F\times G-H)\)。于是答案即为 \(F\times F\times F-3\times (F\times G-H)-H\)。注意最后要除以 \(6\)。
例 4 [BZOJ3160] 万径人踪灭
这道题告诉我们一个事实：FFT 可以实现对字符串的操作。
考虑到回文串的实质实际上就是下标和相等的字符相等。看到下标和相等就可以联想到和卷积，现在的问题是我们怎样构造多项式使得我们能够区分出相等字符和不相等字符。
考虑卷积需要将两项相乘，于是我们令 a 代表 \(1\)，b 代表 \(-1\)，如此积为 \(1\) 代表相等，为 \(-1\) 代表不相等。利用卷积结果可以求出每个回文中心两边有多少对相等字符，设其为 \(cnt\)，则我们直接给答案加上 \(2^{cnt}-1\) 即可。
然后我们现在没有判掉的就是连续的回文串，用 Manacher 再跑一遍即可。
2 快速数论变换 NTT
前置知识：原根相关。
2.1 引入
上文中我们介绍了 FFT，它可以在 \(O(n\log n)\) 复杂度内进行多项式乘法。不过在很多情境下，我们都需要对多项式的结果取模。此时如果继续使用 FFT 的话精度就难以保证了，所以我们需要更换一下 DFT 的方式。
于是我们就引入了数论变换（Number-Theoretic Transform，NTT）来解决模意义下的多项式乘法。
2.2 替换单位根
实际上，在模意义下，FFT 过程中只有单位根不能再使用。我们考虑模意义下有没有能够代替单位根的东西。考虑上面 FFT 中我们运用了单位根的那些性质：
- \(\omega_{n}^{k}=\omega_{2n}^{2k}\)。
- \(\omega_{n}^{k}=\omega_{n}^{k+n}\)。
- \(\omega_{2n}^{n+k}=-\omega_{2n}^k\)。
- \(\omega_n^k\) 对于 \(k\in[0,n-1]\) 互不相等，且 \(\omega_n^n=1\)。
那么什么数可以代替单位根呢？实际上就是上面提到的原根。
对于一个质数 \(p\)，其原根 \(g\) 满足 \(g,g^2,\cdots,g^{p-1}\) 在模 \(p\) 意义下互不相等。为了满足最后一个条件中 \(\omega_n^n=1\) 的条件，我们需要找出一个 \((g^r)^n\equiv 1\pmod p\)，显然只需要让 \(rn=\varphi(p)\) 即可，于是 \(r=\tfrac{\varphi(p)}{n}=\tfrac {p-1}n\)。我们让 \(g^{r}\) 作为新的单位根，来检验一下它是否满足上述性质：
-
\(\omega_n^k=\omega_{2n}^{2k}\)。
\((g^{\frac{p-1}n})^k\equiv g^{\frac{(p-1)k}{n}}\equiv g^{\frac{(p-1)2k}{2n}}\equiv(g^{\frac{p-1}{2n}})^{2k}\pmod p\)。
-
\(\omega_n^k=\omega_n^{k+n}\)。
\((g^{\frac{p-1}{n}})^{n+k}\equiv (g^{\frac{p-1}{n}})^{k}\times g^{p-1}\equiv (g^{\frac{p-1}{n}})^{k}\pmod p\)。
-
\(\omega_{2n}^{n+k}=-\omega_{2n}^k\)。
\((g^{\frac{p-1}{2n}})^{n+k}\equiv (g^{\frac{p-1}{2n}})^{k}\times g^{\frac{p-1}{2}}\pmod p\)。只需证明 \(g^{\frac{p-1}2}\equiv -1\pmod p\) 即可。
由于 \((g^{\frac{p-1}2})^2\equiv g^{p-1}\equiv 1\pmod p\)，于是 \(g^{\frac{p-1}2}\equiv \pm 1\pmod p\)。又由于 \(g^{\frac{p-1}2}\not\equiv g^{p-1}\pmod p\)，所以 \(g^{\frac{p-1}2}\equiv -1\pmod p\)。
-
最后一个条件上文已经说明过，不再赘述。
那么我们就可以确定新的 “单位根” 就是 \(g^{\frac{p-1}n}\)，那么这就要求 \(n\mid p-1\)，由于 \(n\) 就是 \(2^{k}\)，因此我们的 \(p\) 应该是 \(a2^k+1\) 的形式。下面列举几个常见的 NTT 模数及其原根：
- \(998244353=119\times 2^{23}+1,g=3\)。
- \(469762049=7\times 2^{26}+1,g=3\)。
- \(1004535809=479\times 2^{21}+1,g=3\)。
所以我们其实只需要将 FFT 中的 \(\omega_n\) 换成 \(g^r\) 即可，剩下的部分都不用改（对于 IDFT 中求 \(\omega_n^{-k}\) 的部分直接求一个逆元即可）。
2.3 代码
const int Mod = 1004535809;
const int G = 3;//原根
const int InvG = 334845270;//原根逆元
int r[Maxn];
struct Poly {
int n; vector <int> a;
void reset(int len) {n = len, a.resize(len + 1);}
int& operator [](int x) {return a[x];}
void NTT(int len, int typ) {
reset(len - 1);
for(int i = 0; i < len; i++) if(i < r[i]) swap(a[i], a[r[i]]);
for(int h = 1; h < len; h <<= 1) {
int cur = qpow(typ == 1 ? G : InvG, (Mod - 1) / (h << 1), Mod);
for(int i = 0; i < len; i += (h << 1)) {
int w = 1;
for(int j = 0; j < h; j++, w = w * cur % Mod) {
int x = a[i + j], y = w * a[i + j + h] % Mod;
a[i + j] = (x + y) % Mod;
a[i + j + h] = (x - y + Mod) % Mod;
}
}
}
if(typ == -1) {
int iv = qpow(len, Mod - 2, Mod);
for(int i = 0; i < len; i++) a[i] = a[i] * iv % Mod;
}
}
Poly operator * (Poly y) {
Poly x = *this, z;
int n = x.n + y.n, len = 1;
while(len <= n) len <<= 1;
for(int i = 0; i < len; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) * (len >> 1));
x.NTT(len, 1), y.NTT(len, 1);
z.reset(len - 1);
for(int i = 0; i < len; i++) z[i] = x[i] * y[i] % Mod;
z.NTT(len, -1);
z.reset(n);
return z;
}
}F;
2.4 例题
例 1 [SDOI2015] 序列统计
这个题看上去就是一个普通的背包问题，但是权值从和变成了乘。考虑怎样将乘法化为加法，自然可以想到对数。那么在模 \(m\) 意义下的对数自然就是离散对数，也就是指标了。
于是我们对于 \(S\) 中所有数取指标，然后构造多项式 \(F\)，\(a_ix^i\) 项系数是指标为 \(i\) 的数的个数。这样的话我们只需要求出 \(F^n\) 中次数模 \(\varphi(m)=m-1\) 为 \(\gamma(x)\) 的所有项的系数之和即为答案。直接做一遍快速幂即可，同时在乘的时候就可以直接将次数 \(\ge \varphi(m)\) 的部分加到前面然后清空，这样多项式的项数就不会炸了。
例 2 [HAOI2018] 染色
首先满足要求的最多颜色数是 \(\min(M,\tfrac NS)\)。然后我们现在希望求出恰有 \(k\) 种颜色出现次数为 \(S\) 的方案数，不妨记作 \(F(k)\)。那么这个显然是难求的，但是显然我们可以求出至少 \(k\) 种颜色出现次数为 \(S\) 的方案数，记其为 \(G(k)\)，容易得到：
相信这一步读者可以自行完成。然后根据二项式反演可以得到：
将中间的组合数拆开可以得到：
我们发现后面的和式正好可以分为两个部分，并且满足差卷积的形式，因此构造出两个多项式然后直接卷积即可求出 \(F(k)\)，最后算出 \(\sum\limits_{i=0}^{\min(M,\tfrac NS)} F(i)\times W(i)\) 即为答案。
3 任意模数多项式乘法 MTT
上文中介绍的 NTT 只有在模数可以写成 \(a2^k + 1\) 的形式时才能能使用，而如果这个模数没有任何特殊性质，我们就需要辱骂出题人使用任意模数的多项式乘法了。
实现 MTT 的方式有两种：三模数 NTT 或者拆系数 FFT，两种方法的优势区间完全互补。
3.1 三模数 NTT
三模数 NTT，优点是精度大，缺点是常数也大、跑得慢。
顾名思义，三模数 NTT 就是让我们直接计算出多项式的系数取模三个不同的 NTT 模数后的值，然后利用中国剩余定理将其合并；根据中国剩余定理的结论，我们最后合并出来的值应该是取模三个质数乘积的值。而在一般情境下我们的模数是 \(10^9\) 级别，项数是 \(10^5\) 级别，那么我们卷积后得到的最大值应该是 \(10^9\times 10^9\times 10^5=10^{23}\) 级别，而三个模数的乘积可以直接达到 \(10^{27}\) 级别，所以足够我们算出最后的正确结果。
不过这样做的话相当于一次乘法要进行 \(9\) 次 NTT，并且合并的时候还要用 __int128 辅助计算，常数巨大且不美观。所以相比之下这种做法的缺陷要更大一些。
3.2 拆系数 FFT
拆系数 FFT，优点是常数小，缺点是精度可能会炸，要靠 long double 信仰跑。
我们考虑把原多项式的系数按照二进制拆开，把每一个系数写成 \(a\times 2^{15} +b\)。令 \(c=2^{15}\)，则我们可以将 \(f(x),g(x)\) 两个多项式分别写成 \(A(x)+cB(x)\) 和 \(C(x)+cD(x)\) 的形式，此时相乘会发现答案就是：
那么我们只需要求出这四个多项式乘法的值即可。暴力做的话需要进行 \(8\) 次 FFT，效率太低。考虑上文中讲过的 FFT 三次变两次优化，在这里套用相同的思想，我们把 \(C(x),D(x)\) 合并到一个多项式 \(T(x)\) 中，即令 \(T(x)=C(x)+i\times D(x)\)，那么我们只需要求出 \(A(x)T(x)\) 和 \(B(x)T(x)\) 的结果即可，总共做 \(5\) 次 FFT，常数大大降低。
当然你可能会想到把 \(A(x),B(x)\) 也放进同一个多项式，这样就只用做 \(4\) 次 FFT 了。不过这样的话最后求系数会比这个麻烦一点，看个人喜好了。并且实际上只要写法正常，那么拆系数 FFT 也不太会掉精度，实在不行可以预处理单位根。
模板题：【模板】任意模数多项式乘法，采用拆系数 FFT 的代码如下：
#include <bits/stdc++.h>
#define il inline
using namespace std;
typedef long double db;
const int Maxn = (1 << 18) + 5;
const int Inf = 2e9;
const int V = (1 << 15) - 1;
const db Pi = acos(-1.0);
int Mod;
il int Add(int x, int y) {return x + y >= Mod ? x + y - Mod: x + y;} il void pls(int &x, int y) {x = Add(x, y);}
il int Del(int x, int y) {return x - y < 0 ? x - y + Mod : x - y;} il void sub(int &x, int y) {x = Del(x, y);}
il int qpow(int a, int b) {int res = 1; for(; b; a = 1ll * a * a % Mod, b >>= 1) if(b & 1) res = 1ll * res * a % Mod; return res;}
il int Inv(int a) {return qpow(a, Mod - 2);}
template <typename T> il void chkmax(T &x, T y) {x = (x >= y ? x : y);}
template <typename T> il void chkmin(T &x, T y) {x = (x <= y ? x : y);}
template <typename T>
il void read(T &x) {
x = 0; char ch = getchar(); bool flg = 0;
for(; ch < '0' || ch > '9'; ch = getchar()) flg = (ch == '-');
for(; ch >= '0' && ch <= '9'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
flg ? x = -x : 0;
}
template <typename T>
il void write(T x, bool typ = 1) {
static short Stk[50], Top = 0;
x < 0 ? putchar('-'), x = -x : 0;
do Stk[++Top] = x % 10, x /= 10; while(x);
while(Top) putchar(Stk[Top--] | 48);
typ ? putchar('\n') : putchar(' ');
}
il void IOS() {ios::sync_with_stdio(0); cin.tie(0), cout.tie(0);}
il void File() {freopen("in.txt", "r", stdin); freopen("out.txt", "w", stdout);}
bool Beg;
int n, m;
struct comp {
db a, b;
comp operator + (comp y) {return {a + y.a, b + y.b};}
comp operator - (comp y) {return {a - y.a, b - y.b};}
comp operator * (comp y) {return {a * y.a - b * y.b, a * y.b + b * y.a};}
comp operator / (int y) {return {a / y, b / y};}
};
int r[Maxn];
void FFT(comp *a, int len, int typ) {
for(int i = 0; i < len; i++) if(i < r[i]) swap(a[i], a[r[i]]);
for(int h = 1; h < len; h <<= 1) {
comp cur = {cos(Pi / h), typ * sin(Pi / h)};
for(int i = 0; i < len; i += (h << 1)) {
comp w = {1, 0};
for(int j = 0; j < h; j++, w = w * cur) {
comp x = a[i + j], y = a[i + j + h] * w;
a[i + j] = x + y;
a[i + j + h] = x - y;
}
}
}
if(typ == -1) {
for(int i = 0; i < len; i++) a[i] = a[i] / len;
}
}
int real(db x) {return ((long long)(x + 0.5)) % Mod;}
int f[Maxn], g[Maxn];
void MTT(int *f, int *g, int n, int m) {
static comp P[Maxn], Q[Maxn], T[Maxn];
int len = 1; while(len <= n + m) len <<= 1;
for(int i = 0; i <= n; i++) P[i].a = f[i] & V, Q[i].a = f[i] >> 15;
for(int i = 0; i <= m; i++) T[i].a = g[i] & V, T[i].b = g[i] >> 15;
for(int i = 0; i < len; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) * (len >> 1));
FFT(P, len, 1), FFT(Q, len, 1), FFT(T, len, 1);
for(int i = 0; i < len; i++) P[i] = P[i] * T[i], Q[i] = Q[i] * T[i];
FFT(P, len, -1), FFT(Q, len, -1);
for(int i = 0; i < len; i++) f[i] = (real(P[i].a) + (1ll * (real(P[i].b) + real(Q[i].a)) << 15) + (1ll * real(Q[i].b) << 30)) % Mod;
}
bool End;
il void Usd() {cerr << (&Beg - &End) / 1024.0 / 1024.0 << "MB " << (double)clock() * 1000.0 / CLOCKS_PER_SEC << "ms\n"; }
int main() {
read(n), read(m), read(Mod);
for(int i = 0; i <= n; i++) read(f[i]);
for(int i = 0; i <= m; i++) read(g[i]);
MTT(f, g, n, m);
for(int i = 0; i <= n + m; i++) write(f[i], i == n + m);
Usd();
return 0;
}
浙公网安备 33010602011771号