当前位置：首页 > news >正文

多项式全家桶（wjc）

news 2025/8/2 4:53:58

FFT

准备自己写一个高精度的板子，但是不会 FFT 于是学一学 QaQ 。
~~窝还是太菜了，这都不会~~ 。

明确问题

首先，来明确一下我们需要解决的问题：给定一个 $n$ 次多项式 $F(x)$，和一个 $m$ 次多项式 $G(x)$。请求出 $F(x)$ 和 $G(x)$ 的卷积。

前置知识

向量

同时具有大小和方向的量，在几何中通常用带有箭头的线段表示。

圆的弧度制

等于半径长的圆弧所对的圆心角叫做 1 弧度的角，用符号 rad 表示，读作弧度。用弧度作单位来度量角的制度叫做弧度制。

复数

定义

$a,b$ 均为实数，$i$ 为 $\sqrt {-1}$，此时形如 $a+bi$ 的数是复数。其中 $i$ 被称为虚数单位。复数 $a+bi$ 可以由复平面上的由 $(0,0)$ 到 $(a,b)$ 的向量表示。

模长：从原点 $(0,0)$ 至点 $(a,b)$ 的距离，也就是 $\sqrt {a^2+b^2}$ 。
辐角：以逆时针为正方向，$x$ 轴正半轴到给定向量的转角。

共轭复数

$a+bi$ 的共轭复数是 $a-bi$ 。几何意义上，两者模长相同，辐角之和为 $2π$ 。两个共轭复数的和、乘积均为有理数。

复数的加减乘除

$\begin{cases} (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+bdi^2=(ac-bd)+(ad+bc)i\\ \dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)+(bc-ad)i}{c^2+d^2}=(\dfrac{ac+bd}{c^2+d^2})+(\dfrac{bc-ad}{c^2+d^2})i \end{cases}$

几何意义上，复数乘法为模长相乘，辐角相加。

可以轻松写出代码（$x$ 是实部，$y$ 是虚部）：

struct comp{double x, y;comp(double a = 0, double b = 0) {x = a, y = b;}
};
comp operator + (comp x, comp y) {return comp(x.x + y.x, x.y + y.y);}
comp operator - (comp x, comp y) {return comp(x.x - y.x, x.y - y.y);}
comp operator * (comp x, comp y) {return comp(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);}
comp operator / (comp x, comp y) {double a = x.x, b = x.y, c = y.x, d = y.y;return comp((a * c + b * d) / (c * c + d * d), (b * c - a * d) / (c * c + d * d));
}

欧拉公式

$e^{iθ}=\cosθ+i \sin θ$。
在 $θ=\pi$ 时，有 $e^{i\pi}+1=0$。

单位根

定义

$n$ 次单位根指在复数域中，满足 $x^n=1$ 的 $x$。根据代数基本定理，有 $n$ 次单位根共有 $n$ 个，特别的，我们将 $n$ 次本源单位根（这里特指 $e^{i\frac{2\pi}{n}}$）记为 $w_n$。

重要的结论：按辐角大小排列，第 $k$ 个 $n$ 次单位根为 $e^{i\frac{2k\pi}{n}}$ ，也就是 $w_n^k$。

证明：$(e^{i\frac{2k\pi}{n}})^n=e^{i2k\pi}=(e^{i\pi})^{2k}=1$。$e^{i\frac{2k\pi}{n}}$ 是原方程的一个解。
又 $e^{i\frac{2k\pi}{n}}$ 在 $0 \le k < n$ 时互不相同，根据代数基本定理，该方程有且仅有 $n$ 个解。因此该方程解的与 $e^{i\frac{2k\pi}{n}}$ 一一对应。证毕。

引理 & 证明

引理：$w_n^k=w_{nt}^{kt}$ 。

证明：$w_n^k=e^{i \frac{2k\pi}{n}}=e^{i \frac{2kt\pi}{nt}}=w_{nt}^{kt}$ ，证毕。

引理：在 $2\mid n$ 且 $0 \le k < \dfrac{n}{2}$ 时，令 $m=\dfrac{n}{2}$ 有 $w_n^{k+m}=-w^k_n,(w_n^k)^2=(w_n^{k+m})^2=w_m^k$ 。

证明：

$w_n^{k+m}=w_n^k×w_n^m=w_n^k×e^{i\frac{2m\pi}{n}}=w_n^k×e^{i\pi}=-w_n^k$ 。
$(w_n^{k+m})^2=(w_n^k)^2=w_n^{2k}=w_m^k$ 。

证毕。

引理：$\sum_{i=0}^{n-1}w_n^{ik}$ 在 $k=0$ 时值为 $n$ ，否则为 $0$ 。$(0 \le k < n)$

证明：

$k=0$ 时，$w_n^{ik}$ 均为 $1$，故原式值为 $n$ 。
$k≠0$ 时，$w_n^{ik}$ 为公比不为 $1$ 的等比数列。故原式值为 $\dfrac{w^{nk}_ n-w_n^0}{w^k_n-1}=0$ 。

证毕。

解法

FFT & IFFT 能做什么？

FFT（快速傅里叶变换）：在 $O(n\log_2 n)$ 的时间复杂度内，将一个 $n$ 次多项式，从系数表示法，转化为点值表示法。
IFFT（快速傅里叶逆变换）：在 $O(n\log_2 n)$ 的时间复杂度内，将一个 $n$ 次多项式，从点值表示法，转化为系数表示法。

转换为点值表示法后，怎么做乘法？很简单，把每个点值分别乘起来即可。就是酱紫（已经提前定义了复数的乘法）：

for (int i = 0; i < len; i++) Ans[i] = A[i] * B[i];

其中 len 指长度，$A_i$ 和 $B_i$ 指已经算出来的点值，我们把答案存在 $Ans_i$ 里。也可以省去 $Ans$ 数组，写成酱紫：

for (int i = 0; i < len; i++) A[i] = A[i] * B[i];

注：以下均默认 $n$ 为 $2$ 的整数次幂，否则可以在高位填 $0$ 补全。

FFT

递归版 FFT

令 $m=\dfrac{n}{2}$。

$A(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}$。

$\begin{cases} A_0(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^{m-1}\\ A_1(x)=a_1+a_3x+a_5x^2+\dots+a_{n-1}x^{m-1} \end{cases}$

则有 $A(x)=A_0(x^2)+xA_1(x^2)$ 。

在求出 $A_0(x)$ 和 $A_1(x)$ 在 $w_m^0,w_m^1,\dots,w_m^{m-1}$ 的点值后，我们发现，在 $0 \le k < m$ 时：

$\begin{cases} A(w_n^k)=A_0(w_m^k)+w_n^kA_1(w_m^k)\\ A(w_n^{k+m})=A_0(w_m^{k+m})+w_n^{k+m}A_1(w_m^{k+m})=A_0(w_m^k)-w_n^kA_1(w_m^k) \end{cases}$

于是，我们可以通过 $A_0(x)$ 和 $A_1(x)$ 在 $w_m^0,w_m^1,\dots,w_m^{m-1}$ 的点值，$O(n)$ 求出 $A(x)$ 在 $w_n^0,w_n^1,\dots,w_n^{n-1}$ 的点值。这个操作这叫做 蝴蝶操作。

当 $n=1$ 时，$A(w_1^0)=a_0$，直接返回即可。

于是，通过递归，我们就可以 $O(n \log_2n)$ 完成 FFT 了。

迭代版 FFT

递归的常数超级大，过不了，我们需要迭代。怎么实现呢？很简单，既然递归是自上而下的，那么我们的迭代就是一个自下而上的合并过程。我们需要知道递归时调用原多项式系数下标的递归树。

我们其实只需要把递归树的最后一层拿出来，因为我们只要知道最后一层的系数的编号顺序，就可以三重循环：枚举区间长度，枚举区间左端点，枚举位置，完成对区间的合并，算出 $A(x)$ 在 $w_n^0,w_n^1,\dots,w_n^{n-1}$ 的点值了。

这里顺便提一下存储的问题。令 $m=\dfrac{n}{2}$，当我们要求 $A(x)$ 在 $w_n^0,w_n^1,\dots,w_n^{n-1}$ 的点值时， $A_0(x)$ 和 $A_1(x)$ 在 $w_m^0,w_m^1,\dots,w_m^{m-1}$ 的点值已经存储在 $A$ 数组里了。$A_0(w_m^k)$ 存储在 $A_{t+k}$ 里，则 $A_1(w_m^k)$ 存储在 $A_{t+m+k}$ 里（前 $t$ 个位置是已经处理过的）。我们记录下 $A_0(w_m^k)+w_n^kA_1(w_m^k)$ 后，算出 $A(w_n^k)$ 和 $A(w_n^{k+m})$ 把 $A(w_n^k)$ 存储在 $A_{t+k}$ 里，把 $A(w_n^{k+m})$ 存储在 $A_{t+m+k}$ 里。我们发现，因为一次蝴蝶操作刚好是通过原来 $A_{t+k}$ 和 $A_{t+m+k}$ 的值，算出新的 $A_{t+k}$ 和 $A_{t+k+m}$ 的值，可以用类似滚动数组的操作，只用一个数组完成操作。

这一部分的代码：

for (int i = 1; i < Lim; i <<= 1) {comp w = comp(cos(Pi / i), sin(Pi / i));for (int j = 0; j < Lim; j += (i << 1)) {comp w0 = comp(1, 0);for (int k = j; k < i + j; k++, w0 = w * w0) {comp f = A[k], g = w0 * A[k + i];//这也是一个小的常数优化A[k] = f + g, A[k + i] = f - g;}}
}

手玩一下 $n=8$ 的情况，我们发现，递归树的最后一层是：

0 4 2 6 1 5 3 7

把它们转化为二进制：

(000) (100) (010) (110) (001) (101) (011) (111)

把二进制翻转一下：

(000) (001) (010) (011) (100) (101) (110) (111)

再转换为十进制：

0 1 2 3 4 5 6 7

诶，我们发现递归树的最后一层的第 $i$ 个数 $num_i$ 其实就是把 $i$ 的二进制翻转一下的结果。至于为什么是这样子的……本蒟蒻不会证。但感性理解一下，一层层递归大概相当于往二进制里从右到左填数，递归到 $A_0$ 一侧的就在那一位填 $0$，递归到 $A_1$ 那一侧的就在那一位填 $1$ 。当然如果直接反转常数会爆表，所以，我们需要其他方式求出 $num_i$。有下面一个递推式 $num_i = \left\lfloor\dfrac{num_{\left\lfloor\frac{i}{2}\right\rfloor}}{2}\right\rfloor+[i\bmod 2=1]\times\dfrac{n}{2}$。这个公式不难自证，我就不在此赘述了。这里注意，我们虽然在函数的开始对 $A$ 数组进行了交换，但由于上述推导，函数结束后 $A$ 数组为按 $w_n ^ 0, w_n ^ 1, \cdots, w_n ^ {n - 1}$ 顺序的 $A(x)$ 点值。

FFT 完整代码：

inline void FFT(int Lim, comp *A) {for (int i = 0; i < Lim; i++)if (i < num[i]) swap(A[i], A[num[i]]);for (int i = 1; i < Lim; i <<= 1) {comp w = comp(cos(Pi / i), sin(Pi / i));for (int j = 0; j < Lim; j += (i << 1)) {comp w0 = comp(1, 0);for (int k = j; k < i + j; k++, w0 = w * w0) {comp f = A[k], g = w0 * A[k + i];A[k] = f + g, A[k + i] = f - g;}}} 
}

IFFT

已知 $b_k=A(w_n^k)=\sum_{i=0}^{n-1}a_i(w_n^k)^i$，求 $a_k$ 。

\[\begin{aligned} B(w_n ^ {n - k}) & = \sum_{i = 0} ^ {n - 1}b_i(w_n ^ {n - k}) ^ i \\ & = \sum_{i = 0} ^ {n - 1} \left(\sum_{j = 0} ^ {n - 1}a_j(w_n ^ i) ^ j \right)(w_n ^ {-k}) ^ i \\ & = \sum_{i = 0} ^ {n - 1} \sum_{j = 0} ^ {n - 1}a_j(w_n ^ i) ^ j(w_n ^ {-k}) ^ i\\ & = \sum_{j = 0} ^ {n - 1}a_j\left(\sum_{i = 0} ^ {n - 1}(w_n ^ i) ^ j(w_n ^ {-k}) ^ i\right)\\ & = \sum_{j = 0} ^ {n - 1}a_j\left(\sum_{i = 0} ^ {n - 1}w_n ^ {ij - ik}\right)\\ & = \sum_{j = 0} ^ {n - 1}a_j\left(\sum_{i = 0}^{n - 1}w_n ^ {i(j - k)}\right) \end{aligned} \]

$j=k$ 时，$\sum_{i=0}^{n-1}w_n^{i(j-k)}=n$ ，否则，由于 $|j-k|<n$ ，$\sum_{i=0}^{n-1}w_n^{i(j-k)}=0$，所以原式的值为 $na_k$ 。

$a_k=\dfrac{1}{n}B(w_n^{n-k})$ ，也就是说，我们只要求出 $B(x)$ 在 $w_n^0,w_n^1,\dots,w_n^{n-1}$ 的点值，然后将这些点值全部除以 $n$ ，最后将 $a_1$ ~ $a_{n-1} $ 翻转即可（注意是 $a_1$ 和 $a_{n-1}$ 交换，而不是 $a_0$ 和 $a_{n-1}$ 交换）。这可以在一次 FFT 后轻松 $O(n)$ 解决。于是，IFFT 结束！上代码。

inline void IFFT(int Lim, comp *A) {FFT(Lim, A); reverse(A + 1, A + Lim);for (int i = 0; i < Lim; i++) A[i].x /= Lim, A[i].y /= Lim;
}

IFFT 还有一种更常见的写法，就是把 FFT 和 IFFT 合成一个函数，用一个 flag 表示是 FFT 还是 IFFT。flag 为 $1$ 时就是 FFT，为 $-1$ 时是 IFFT。上代码。

inline void FFT(int Lim, comp *A, int flag) {for (int i = 0; i < Lim; i++)if (i < num[i]) swap(A[i], A[num[i]]);for (int i = 1; i < Lim; i <<= 1) {comp w = comp(cos(Pi / i), sin(Pi / i) * flag);for (int j = 0; j < Lim; j += (i << 1)) {comp w0 = comp(1, 0);for (int k = j; k < i + j; k++, w0 = w * w0) {comp f = A[k], g = w0 * A[k + i];A[k] = f + g, A[k + i] = f - g;}}} 
}

这个执行完 IFFT 之后，输出之前，还要除以 $n$，这一步没有写在函数里。

三次变两次优化

由于 $A(x)$ 和 $B(x)$ 均为整系数多项式，我们可以把 $B(x)$ 塞到 $A(x)$ 的虚部上去，就是令 $F(x)=A(x)+iB(x)$，然后算出 $F(x)^2$，这共需要两次 FFT，$F(x)^2=(A(x)+iB(x))^2=(A(x)^2-B(x)^2)+i(2A(x)B(x))$。所以，我们只要把 $F(x)^2$ 的虚部取出来，再除以二，就能算出 $A(x)$ 和 $B(x)$ 的卷积了。

模板题代码：

#include <bits/stdc++.h>
#define LL long long
using namespace std;namespace Polynomial{const double Pi = acos(-1), eps = 0.5;const int MAXN = 4000005, p = 998244353, g = 3, invg = 332748118;struct comp{double x, y;comp(double a = 0, double b = 0) {x = a, y = b;}};comp operator + (comp x, comp y) {return comp(x.x + y.x, x.y + y.y);}comp operator - (comp x, comp y) {return comp(x.x - y.x, x.y - y.y);}comp operator * (comp x, comp y) {return comp(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);}comp operator / (comp x, comp y) {double a = x.x, b = x.y, c = y.x, d = y.y;return comp((a * c + b * d) / (c * c + d * d), (b * c - a * d) / (c * c + d * d));}int num[MAXN];comp A[MAXN], B[MAXN], C[MAXN];inline void FFT(int Lim, comp *A, int flag) {for (int i = 0; i < Lim; i++)if (i < num[i]) swap(A[i], A[num[i]]);for (int i = 1; i < Lim; i <<= 1) {comp w = comp(cos(Pi / i), sin(Pi / i) * flag);for (int j = 0; j < Lim; j += (i << 1)) {comp w0 = comp(1, 0);for (int k = 0; k < i; k++, w0 = w * w0) {comp f = A[j + k], g = A[i + j + k] * w0;A[j + k] = f + g, A[i + j + k] = f - g;}}}if (flag == -1)for (int i = 0; i < Lim; i++) A[i] = A[i] / Lim;}inline void Multiply(LL *a, LL *b, int n, int m) {for (int i = 0; i <= n; i++) A[i].x = a[i];for (int i = 0; i <= m; i++) A[i].y = b[i];int len = 1; while (len <= n + m) len <<= 1;for (int i = 0; i < len; i++) num[i] = (num[i >> 1] >> 1) | (i & 1) * (len >> 1);FFT(len, A, 1); for (int i = 0; i < len; i++) A[i] = A[i] * A[i]; FFT(len, A, -1);for (int i = 0; i <= n + m; i++) a[i] = (LL)(A[i].y / 2 + eps);}
}
using namespace Polynomial; int n, m;
LL a[MAXN], b[MAXN];
signed main() {cin >> n >> m;for (int i = 0; i <= n; i++) cin >> a[i];for (int i = 0; i <= m; i++) cin >> b[i];Multiply(a, b, n, m);for (int i = 0; i <= n + m; i++) cout << a[i] << " ";return 0;
}

NTT

众所周知，FFT 常数大，精度低，所以说，我们可能需要一些更高科技的东西，也就是 NTT！

首先明确很重要的一点，FFT 里的 $w_i$ 可以不指单位根。我们要实现 FFT，仅需要 $w_i$ 这种符号满足上文中所提到的，单位根满足的三个性质。

阶

若 $\gcd(a,p)=1$ 且 $p>1$，对于满足 $a^x\equiv1\pmod{p}$ 的最小正整数 $x$，称为 $a$ 模 $p$ 的阶，记作 $\delta_p(a)$。由欧拉定理，在 $\gcd(a,p)=1$ 时，有 $a^{\varphi(p)}\equiv1\pmod{p}$，所以，在 $\gcd(a,p)=1$ 时，$a$ 模 $p$ 的阶存在。

性质：$a,a^2,\dots,a^{\delta_p(a)}$ 模 $p$ 两两不同余。

证明：若有 $i,j$ 满足 $1\le i<j\le \delta_p(a)$ 且 $a^i\equiv a^j\pmod{p}$，则有 $a^{j-i}\equiv1\pmod{p}$，又 $1\le j-i<\delta_p(a)$，这与阶的定义矛盾，故：$a,a^2,\dots,a^{\delta_p(a)}$ 模 $p$ 两两不同余。

当然，阶还有很多其他的性质，在此略去，可以在 oi-wiki 上查看详细的性质、证明。

原根

设 $p$ 为正整数，$a$ 为整数，若 $\gcd(a,p)=1$ 且 $\delta_p(a)=\varphi(p)$，则称 $a$ 为模 $p$ 的原根。

栗子：由于 $\delta_{10}(3)=4=\varphi(10)$，所以 $3$ 是模 $10$ 的原根。这个例子也证明了不仅仅是质数有原根。

有以下两个结论：

如果一个数 $m$ 有原根，则它的原根个数是 $\varphi(\varphi(m))$。
$m$ 有原根的充分必要条件为 $m=2,4,p^\alpha,2\times p^\alpha$，其中 $p$ 为奇素数，$\alpha$ 为正整数。

这两条结论的证明以及其他结论见 oi-wiki。

NTT中的 $w_i$

直接放结论，在 NTT 中，$w_i=g^{\frac{p-1}{i}}$，其中 $p$ 是模数（质数），$g$ 是模 $p$ 的原根。让我们来验证一下此时的 $w_i$ 是否满足上文中提及的三个性质。

引理：$w_n^k=w_{nt}^{kt}$ 。

证明：$w_n^k=g^{\frac{p-1}{n}k}=g^{\frac{p-1}{nt}kt}=w_{nt}^{kt}$ ，证毕。

引理：在 $2\mid n$ 且 $0 \le k < \dfrac{n}{2}$ 时，令 $m=\dfrac{n}{2}$ 有 $w_n^{k+m}\equiv-w^k_n\pmod{p},(w_n^k)^2\equiv(w_n^{k+m})^2\equiv w_m^k\pmod{p}$ 。

证明：
$w_n^n\equiv g^{p-1}\equiv1\pmod{p}$，所以 $w_n^m\equiv \pm1\pmod{p}$。又根据原根的定义， $w_n^m\equiv g^{\frac{p-1}{2}}\not\equiv 1\pmod{p}$，得 $w_n^m\equiv -1\pmod{p}$。所以有：

$w_n^{k+m}\equiv w_n^m\times w_n^k\equiv-w_n^k\pmod{p}$。
$(w_n^{k+m})^2\equiv (w_n^k)^2\equiv w_n^{2k}\equiv w_m^k\pmod{p}$ 。

证毕。

引理：$\left(\sum_{i=0}^{n-1}w_n^{ik}\right)\bmod p$，$k=0$ 时值为 $n$ ，否则为 $0$ 。$(0 \le k < n)$

证明：

$k=0$ 时，$w_n^{ik}$ 模 $p$ 均为 $1$，故原式值为 $n$ 。
$k≠0$ 时，$w_n^{ik}$ 为公比不为 $1$ 的等比数列。故原式值为 $\dfrac{w^{nk}_ n-w_n^0}{w^k_n-1}\equiv 0\pmod{p}$ 。

证毕。

所以，$w_i=g^{\frac{p-1}{i}}$ 时，$w_i$ 满足全部三个性质。所以，我们就可以解决 NTT 了吗？不，还有最后一点：** $\dfrac{p-1}{i}$ 不一定是整数**。

为了使计算中 $\frac{p-1}{i}$ 是整数，我们选定的模数必须为 $a\times2^k+1$，其中 $k \ge 21$ （选 $21$ 是因为 $2^{21} \ge 2\times 10^6$），而 $998244352=119\times2^{23}$，所以，在 OI 里，我们常常使用 $998244353$ 做 NTT 的模数，它有原根 $3$。

最后的最后，还要注意一点，在 IFFT 的最后，我们需要除以 $n$，在 NTT 里，就变成了乘 $n$ 模 $p$ 的逆元。

NTT 结束，上代码：

#include <bits/stdc++.h>
#define int long long
using namespace std;const int MAXN = 3000005;
const int p = 998244353, g = 3, invg = 332748118;
int n, m, len = 1, a[MAXN], b[MAXN], num[MAXN];
inline int ksm(int x, int y, int mod) {int ans = 1;while (y) {if (y & 1) ans = (ans * x) % mod;x = (x * x) % mod;y >>= 1;}return ans % mod;
}inline void NTT(int *A, int Lim, int type) {for (int i = 0; i < Lim; i++)num[i] = (num[i >> 1] >> 1) | (i & 1) * (Lim >> 1);for (int i = 0; i < Lim; i++)if (i < num[i]) swap(A[i], A[num[i]]);for (int i = 1; i < Lim; i <<= 1) {int w = ksm((type == 1 ? g : invg), (p - 1) / (i << 1), p);for (int j = 0; j < Lim; j += (i << 1)) {int w0 = 1;for (int k = j; k < j + i; k++, w0 = (w * w0) % p) {int f = A[k], g = (A[i + k] * w0) % p;A[k] = (f + g) % p, A[i + k] = (f - g + p) % p;}}}if (type == -1) {int inv = ksm(Lim, p - 2, p);for (int i = 0; i < Lim; i++) A[i] = (A[i] * inv) % p;}
}signed main() {cin >> n >> m;for (int i = 0; i <= n; i++) cin >> a[i];for (int i = 0; i <= m; i++) cin >> b[i];while (len <= n + m) len <<= 1;NTT(a, len, 1), NTT(b, len, 1);for (int i = 0; i < len; i++) a[i] = (a[i] * b[i]) % p;NTT(a, len, -1);for (int i = 0; i <= n + m; i++) cout << a[i] << " ";return 0;
}

任意模数多项式乘法

由于 NTT 只能处理形如 $a \times 2 ^ k + 1$ 且 $k$ 较大的模数 $p$，且 FFT 会掉精度，我们需要特殊的方式处理任意模数多项式乘法。我们设 $A(x) = C(x)t + D(x), B(x) = E(x)t + F(x)$，其中 $t = 32768$。我们考虑 $(C(x) + iD(x))E(x)$ 与 $(C(x) + iD(x))F(x)$，我们发现分离两者的实部与虚部可得到 $C(x)E(x), D(x)E(x), C(x)F(x), D(x)F(x)$，从而算出 $(C(x)t + D(x))(E(x)t + F(x))$。为此，我们只需对 $C(x) + iD(x), E(x), F(x)$ 做 $3$ 次 FFT，对于 $(C(x) + iD(x))E(x), (C(x) + iD(x))F(x)$ 做两次 IFFT 即可。代码如下：

#include <bits/stdc++.h>
#define LL long long
using namespace std;namespace Polynomial{const double Pi = acos(-1), eps = 0.5;const int MAXN = 4000005, p = 998244353, g = 3, invg = 332748118;struct comp{double x, y;comp(double a = 0, double b = 0) {x = a, y = b;}};comp operator + (comp x, comp y) {return comp(x.x + y.x, x.y + y.y);}comp operator - (comp x, comp y) {return comp(x.x - y.x, x.y - y.y);}comp operator * (comp x, comp y) {return comp(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);}comp operator / (comp x, comp y) {double a = x.x, b = x.y, c = y.x, d = y.y;return comp((a * c + b * d) / (c * c + d * d), (b * c - a * d) / (c * c + d * d));}int num[MAXN];comp A[MAXN], B[MAXN], C[MAXN], Wn[MAXN];inline void FFT(int Lim, comp *A, int f) {for (int i = 0; i < Lim; i++)if (i < num[i]) swap(A[i], A[num[i]]);for (int i = 1; i < Lim; i <<= 1) for (int j = 0; j < Lim; j += (i << 1)) for (int k = j; k < i + j; k++) {comp w0 = Wn[Lim / i * (k - j)];w0.y *= f;comp f = A[k], g = A[i + k] * w0;A[k] = f + g, A[i + k] = f - g;}}inline void Mul(LL *a, LL *b, int n, int m) {for (int i = 0; i <= n; i++) A[i].x = a[i];for (int i = 0; i <= m; i++) A[i].y = b[i];int len = 1; while (len <= n + m) len <<= 1;for (int i = 0; i < len; i++) Wn[i] = comp(cos(Pi * i / len), sin(Pi * i / len));for (int i = 0; i < len; i++) num[i] = (num[i >> 1] >> 1) | (i & 1) * (len >> 1);FFT(len, A, 1); for (int i = 0; i < len; i++) A[i] = A[i] * A[i]; FFT(len, A, -1);for (int i = 0; i <= n + m; i++) a[i] = (LL)(A[i].y / (2 * len) + eps);}inline void MTT(LL *a, LL *b, int n, int m, int mod) {const int t = 32768;for (int i = 0; i <= n; i++) A[i] = comp(a[i] / t, a[i] % t);for (int i = 0; i <= m; i++) B[i] = b[i] / t, C[i] = b[i] % t;int len = 1; while (len <= n + m) len <<= 1;for (int i = 0; i < len; i++) Wn[i] = comp(cos(Pi * i / len), sin(Pi * i / len));for (int i = 0; i < len; i++) num[i] = (num[i >> 1] >> 1) | (i & 1) * (len >> 1);FFT(len, A, 1), FFT(len, B, 1), FFT(len, C, 1);for (int i = 0; i < len; i++) B[i] = A[i] * B[i], C[i] = A[i] * C[i];FFT(len, B, -1), FFT(len, C, -1);for (int i = 0; i < len; i++) {a[i] = ((LL)(B[i].x / len + eps) % mod * (t * t % mod)) % mod;a[i] = (a[i] + (((LL)(B[i].y / len + eps) + (LL)(C[i].x / len + eps)) % mod * t) % mod) % mod;a[i] = ((LL)(C[i].y / len + eps) % mod + a[i]) % mod;}}
}
using namespace Polynomial; int n, m, k;
LL a[MAXN], b[MAXN];
signed main() {cin >> n >> m >> k;for (int i = 0; i <= n; i++) cin >> a[i];for (int i = 0; i <= m; i++) cin >> b[i];MTT(a, b, n, m, k);for (int i = 0; i <= n + m; i++) cout << a[i] << " ";return 0;
}

分治FFT

P4721

考虑 CDQ 分治，对于区间 $[l, r]$，我们只需考虑 $[l, mid]$ 对 $[mid + 1, r]$ 的贡献。容易发现，贡献为 $f_l, f_{l + 1}, \cdots, f_{mid}$ 与 $g_1, g_2, \cdots, g_{r - l}$ 的卷积后对应项，直接卷积即可。复杂度 $\mathcal{O}(n \log ^ 2 n)$。

UOJ041

题意简述：已知 $f_1 = 1, f_i = f_{i - 1} + \sum\limits_{j = 1} ^ {i - 1} (j - 1) f_j f_{i - j}$，求 $f_1, f_2, \cdots, f_n$。

考虑 $\max(i, j) \in [L, mid]$ 的二元组 $(i, j)$ 对 $i + j$ 的贡献。

$L = 1$ 时，有 $i, j \in [L, mid]$ 将 $f_1, f_2, \cdots, f_{mid}$ 与 $0f_1, 1f_2, \cdots, (mid - 1)f_{mid}$ 卷积即可。
$L \neq 1$ 时，我们发现 $R - L < L$，故 $i, j$ 仅有一个在 $[L, mid]$ 中，将 $f_L, \cdots, f_{mid}$ 与 $0f_1, 1f_2, \cdots, (R - L - 1) f_{R - L}$ 卷积，再将 $(L - 1)f_L, Lf_{L + 1}, \cdots, (mid - 1)f_{mid}$ 与 $f_1, f_2, \cdots, f_{R - L}$ 卷积，再在对应位置贡献即可。

故分治 FFT 即可，复杂度为 $\mathcal{O}(n \log ^ 2 n)$。

全家桶

多项式求逆

明确问题

给定多项式 $A(x)$，求多项式 $B(x)$，满足 $A(x)B(x) \equiv 1 \pmod {x ^ n}$。

做法

我们设 $B_k(x)$ 为 $B(x) \bmod x ^ k$ 的值。考虑倍增，起点为 $B_1(x) = a_0 ^ {-1}$。

假设我们已经求得 $A(x)B_k(x) \equiv 1 \pmod{x ^ k}$，

\[\begin{aligned} & A(x)B_k(x) - 1 \equiv 0 \pmod{x ^ k} \\[1pt] & \implies (A(x)B_k(x) - 1) ^ 2 \equiv 0 \pmod{x ^ {2k}} \\[1pt] & \implies - A(x) ^ 2B_k(x) ^ 2 + 2A(x)B_k(x) \equiv 1 \pmod{x ^ {2k}} \\[1pt] & \implies A(x)(2B_k(x) - A(x)B_k(x) ^ 2) \equiv 1 \pmod{x ^ {2k}} \\[1pt] & \implies B_{2k}(x) \equiv 2B_k(x) - A(x)B_k(x) ^ 2 \pmod{x ^ {2k}} \end{aligned} \]

用 FFT 进行多项式乘法即可，复杂度为 $\mathcal{O}(n \log n)$。

多项式牛顿迭代

前置知识：泰勒展开

$f(k) = \sum\limits_{i \ge 0} \dfrac{f ^ {(i)}(k)}{i!}(x - k) ^ i$，其中 $f ^ {(i)}$ 为 $f$ 的 $i$ 阶导数。

理解方式（不严谨）：设 $f(x) = \sum\limits_{i \ge 0} a_i(x - k) ^ i$，有 $f ^ {(i)} (x) = \sum\limits_{j \ge i} a_j j ^ {\underline{i}}(x - k) ^ {j - i}$，故 $f ^ {(i)}(k) = a_i i!$，带入即得。

明确问题

对于函数 $F(G)$，其中 $G$ 是多项式，我们想求出 $F$ 的零点。

做法

对 $F(G)$ 进行泰勒展开，我们有 $F(G) = F(H) + (G - H) F'(H) + (G - H) ^ 2 F''(H) + \cdots$，其中 $H$ 为任意多项式。我们设 $A$ 为 $F$ 的零点，仍然考虑倍增，假设我们已经求出 $A_n$（即 $A(x) \bmod x ^ n$），考虑求 $A_{2n}$。由泰勒展开，我们有 $F(A) = F(A_n) + (A - A_n)F'(A_n) + (A - A_n) ^ 2F''(A_n) + \cdots$。

\[\begin{aligned} & A_n \equiv A \pmod{x ^ n} \\[1pt] & \implies A_n - A \equiv 0 \pmod {x ^ n} \\[1pt] & \implies (A_n - A) ^ k \equiv 0 \pmod {x ^ {2n}} \ \ \ (k \ge 2) \\[1pt] & \implies 0 = F(A) \equiv F(A_n) + (A - A_n) F'(A_n) \pmod {x ^ {2n}} \\[1pt] & \implies -F(A_n) \equiv (A_{2n} - A_n)F'(A_n) \pmod {x ^ {2n}} \\[1pt] & \implies A_{2n} \equiv A_n - \dfrac{F(A_n)}{F'(A_n)} \pmod {x ^ {2n}} \end{aligned} \]

使用多项式求逆与多项式乘法即可。

多项式开根

前置知识：Cipolla 算法

更严谨的证明见此处。

Cipolla 算法可以对 $\bmod p$ 意义下的二次剩余 $n$，求出 $x$ 满足 $x ^ 2 \equiv n \pmod p$。

考虑随机 $a$ 满足 $a ^ 2 - n$ 不是 $\bmod p$ 意义下的二次剩余，根据二次剩余相关知识，我们知道对于素数 $p$，$x$ 是 $\bmod p$ 意义下的二次剩余 $\iff x ^ {\frac{p - 1}{2}} \equiv 1 \pmod p$，且 $\bmod p$ 意义下共有 $\dfrac{p - 1}{2}$ 个二次剩余，$\dfrac{p - 1}{2}$ 个二次非剩余，所以我们期望 $\mathcal{O}(1)$ 次即可寻找到满足条件的 $a$，且可以 $\mathcal{O}(\log p)$ 判断 $a$ 是否合法。

考虑扩域，定义 $\mathbf{F}_{p ^ 2} = \mathbf{F}_p(\sqrt {a ^ 2 - n})$，即 $\mathbf{F_{p ^ 2}}$ 中的每个数都可以写成 $i + j \sqrt {a ^ 2 - n} \ \ \ (i,j \in \{0, 1, 2, \cdots p - 1\})$ 的形式。对于 $\mathbf{F}_{p ^ 2}$ 中数的运算，可以用类似复数运算的方式解决。下证明 $x \equiv (a + \sqrt{a ^ 2 - n}) ^ {\frac{p + 1}{2}} \pmod p$。令 $w = \sqrt{a ^ 2 - n}$，先证两个性质：

$w ^ p \equiv -w \pmod{p}$：$w ^ p \equiv (a ^ 2 - n) ^ {\frac{p - 1}{2}}w \equiv -w \pmod{p}$。
$(a + b) ^ p \equiv a ^ p + b ^ p \pmod{p}$：注意到对 $1 \le k \le p - 1$，有 $\dbinom{p}{k} \equiv 0 \pmod {p}$，二项式定理展开即证。

只需证 $\Big((a + \sqrt {a ^ 2 - n}) ^ {\frac{p + 1}{2}} \Big) ^ 2 \equiv n \pmod p$，这是因为 $(a + w) ^ {p + 1} \equiv (a ^ p + w ^ p)(a + w) \equiv (a - w)(a + w) \equiv a ^ 2 - w ^ 2 \equiv n ^ 2 \pmod p$，当然，这个证明并不是很严谨，严谨的证明请参考以上链接。