多项式基础 - 学习笔记

1. FFT

前置知识：复数

基础定义

令 \(i\) 为方程 \(x^2 = -1\) 的解

我们定义，形如 \(z = \bm{a+bi}\)，其中 \(a, b \in \mathbb{R}\) 的数 \(z\)，称为复数；对于 \(z = a+bi\)，我们称 \(a\) 为其实部，\(bi\) 为其虚部

定义复数 \(z_1 = z_2\)，当且仅当 \(\begin{cases} a_1 = a_2 \\ b_1 = b_2 \end{cases}\)；由此：

复数 \(z = a+bi\) 唯一对应平面直角坐标系中一点 \((a, b)\)；换句话说，复数集与平面直角坐标系中的点一一对应
复数 \(z = a+bi\) 也唯一对应向量 \(\overrightarrow{OZ}\)，其中点 \(Z\) 的坐标为 \((a, b)\)

我们称对应复数集的平面直角坐标系为复平面，x 轴为实轴，y 轴为虚轴

类似向量，定义复数 \(z = a+bi\) 的模为 \(\bm{|z| = \sqrt{a^2+b^2}}\)

四则运算

令 \(z_1 = a_1+b_1i\)，\(z_2 = a_2+b_2i\)，定义：

复数加法 \(z_1+z_2 = (a_1+a_2)+(b_1+b_2)i\)
复数减法 \(z_1-z_2 = (a_1-a_2)+(b_1-b_2)i\)
复数乘法 \(z_1z_2 = (a_1+b_1i)(a_2+b_2i) = a_1a_2+b_1b_2i^2+(a_1b_2+b_1a_2)i\)
带入 \(i^2 = -1\)，可得 \(z_1z_2 = (a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i\)
复数除法 \(\displaystyle \frac{z_1}{z_2} = \frac{a_1+b_1i}{a_2+b_2i}\)
为使分母化为实数，上下同乘 \(a_2-b_2i\)，可得 \(\displaystyle \frac{z_1}{z_2} = \frac{(a_1+b_1i)(a_2-b_2i)}{(a_2+b_2i)(a_2-b_2i)} = \frac{a_1a_2+b_1b_2}{a_2^2+b_2^2} + \frac{a_2b_1-a_1b_2}{a_2^2+b_2^2}i\)

手写结构体代码如下：

#define db double
struct cp{
	db x, y; cp(db xx = 0, db yy = 0){ x = xx; y = yy; }
	cp operator + (const cp& b){ return cp(x+b.x, y+b.y); }
	cp operator - (const cp& b){ return cp(x-b.x, y-b.y); }
	cp operator * (const cp& b){ return cp(x*b.x-y*b.y, x*b.y+y*b.x); }
	cp operator / (const cp& b){ db t = (b.x*b.x+b.y*b.y); return cp((x*b.x+y*b.y)/t, (y*b.x-x*b.y)/t); }
}

辐角与辐角主值

非零复数 \(z = a+bi\) 对应点 \((a, b)\)，我们用极坐标将其表示为 \((r, \theta)\)，其中 \(\bm{r = |z|}\)

记向量与实轴的夹角 \(\theta\) 为 \(z\) 的辐角，记为 \(\theta = \text{Arg}\ z\)，显然有 \(\displaystyle \tan \theta = \frac{b}{a}\)

注意到任一非零复数 \(z\) 都有无穷多个辐角，我们用 \(\text{arg}\ z\) 表示其中满足 \(-\pi < \theta < \pi\) 的特定 \(\theta\)，称为辐角主值

由此，可以引出复数相乘的重要规律：复数相乘，模长相乘，辐角相加 (非常重要)

复数相乘演示图

证明：设 \(z_1 = a_1+b_1i\)，\(z_2 = a_2+b_2i\)：

先证模长相乘，\(\text{LHS} = \sqrt{a_1^2+b_1^2} \times \sqrt{a_2^2+b_2^2}\)，\(\text{RHS} = \sqrt{(a_1a_2-b_1b_2)^2+(a_1b_2+a_2b_1)^2}\)，展开，显然有 \(\text{LHS} = \text{RHS}\)

再证辐角相加，设 \(z_1 = \overrightarrow{OA}\)，\(z_2 = \overrightarrow{OB}\)，相乘结果为 \(\overrightarrow{OC}\)；取 \(P(1, 0)\)，考虑证 \(\triangle OAP\) 与 \(\triangle OCB\) 相似

由模长相乘，有 \(\overline{OA} \times \overline{OB} = \overline{OC}\)；又显然有 \(\overline{OP} \times \overline{OB} = \overline{OB}\)，那么只需证 \(\overline{AP} \times \overline{OB} = \overline{CB}\)

套下两点距离公式暴力展开即可，不难证得相似

那么有 \(\angle COP = \angle BOP + \angle BOC = \angle BOP + \angle AOP\)，得证

单位根

我们定义满足 \(\bm{|z| = 1}\) 的复数 \(z\) 构成单位圆

称方程 \(x^n = 1\) 在复数意义下的解为 \(n\) 次单位根；这样的解必恰有 \(n\) 个，且 \(n\) 次单位根将单位圆 \(n\) 等分

证明：

由模长相乘，\(|z| < 1\) 时会越乘越小，\(|z| > 1\) 时会越乘越大；于是单位根必然在单位圆上

由辐角相加，有 \(n \theta = 2k \pi\ (k \in \mathbb{Z})\)

取 \(\displaystyle \theta = \frac{1(2 \pi)}{n}, \frac{2(2 \pi)}{n}, \frac{3(2 \pi)}{n}, \cdots, \frac{(n-1)(2 \pi)}{n}\)，至少有这 \(n\) 个解

由代数基本定理，有且仅有这 \(n\) 个解；显然它们将单位圆 \(n\) 等分，证毕

我们记这些单位根为 \(\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}\) (以下均省略指数上的 \(\bmod n\))，它们有如下性质：

\(\omega_{n}^{0} = \omega_{n}^{n} = 1\)
\(\omega_{n}^{k} = (\omega_{n}^{1})^k\)，由模长相乘辐角相加，显然
\(\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i+j}\)，先转 \(i\) 次再转 \(j\) 次即为转 \(i+j\) 次
\(\bm{\omega_{n}^{k} = \omega_{2n}^{2k}}\)，即将每段细分为两份 (重要)
\(\bm{\omega_{2n}^{k+n} = -\omega_{2n}^{k}}\)，由 \(\omega_{2n}^{k+n} = \omega_{2n}^{k} \times \omega_{2n}^{n}\)，多乘 \(\omega_{2n}^{n}\) 相当于多转半圈，关于原点对称 (重要)

DFT 与 IDFT 思想

P.S. 下面使用 \(F_{i}\) 表示多项式 \(F(x)\) 的 \(i\) 次项系数

两多项式相乘(卷积)定义为 \(H(x) = \displaystyle \sum_{i=0} \sum_{j=0} F_iG_jx^{i+j}\)

令 \(n\) 次多项式 \(F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\)；这称为多项式的系数表达

系数表达较为常用；但系数表达直接算卷积是 \(O(n^2)\) 的
由 \(n\) 个等式可以解出 \(n\) 元方程的解，我们可以用 \(n+1\) 个点唯一确定 \(F(x)\)；这称为多项式的点值表达

这种表达方法的好处是：
- 若我们知道 \(F(x)\) 与 \(G(x)\) 在 \(x_0\) 处的点值 \(F(x_0)\) 与 \(G(x_0)\)，直接就有 \(H(x_0) = F(x_0) \times G(x_0)\)
  
  换句话说，若我们知道 \(F(x)\) 与 \(G(x)\) 的点值表达，可以 \(O(n)\) 算出 \(H(x)\) 的点值表达
- 系数表达允许我们带入具有特殊性质的值，这有助于创造好的性质

由上，考虑使用点值表达解决问题

那么，现在只缺将系数表达转化为点值表达(求值)，以及将点值表达转化为系数表达(插值)的算法

朴素系数转点值算法即为 DFT，朴素点值转系数算法即为 IDFT

FFT 原理

将朴素 DFT 算法优化，可得 FFT

FFT 的主要思想是利用分治来优化时间复杂度

设有 \(n = 2^k\ (k \in \mathbb{N^*})\)，对于项数不足 \(2^k\) (次数不足 \(2^k-1\) ) 的多项式，我们在其高次项补零即可

将多项式 \(F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\) 根据次数奇偶拆成两半，分别为：

\(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2}\)
\(a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1}\)

我们根据这两部分构造两个多项式 \(FL(x)\) 与 \(FR(x)\)，具体的：

\(FL(x) = a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}\)
\(FR(x) = a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1}\)

那么显然有关系式：\(\displaystyle \bm{F(x) = FL(x^2)+xFR(x^2)}\)

考虑带入 \(\omega_{n}^{k}\)，发掘性质：

\(\displaystyle F(\omega_{n}^{k}) = FL((\omega_{n}^{k})^2)+\omega_{n}^{k}FR((\omega_{n}^{k})^2)\)

\(\displaystyle = FL(\omega_{n}^{2k})+\omega_{n}^{k}FR(\omega_{n}^{2k})\)

\(\displaystyle = FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})\)

由此，只要知道 \(FL(x)\) 与 \(FR(x)\) 在 \(\displaystyle \omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) 处的取值，就可以得到 \(F(x)\) 在 \(\displaystyle \omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{\frac{n}{2}-1}\) 处的取值

对于后半段，考虑带入 \(\displaystyle \omega_{n}^{k+\frac{n}{2}}\)：

\(\displaystyle F(\omega_{n}^{k+\frac{n}{2}}) = FL((\omega_{n}^{k+\frac{n}{2}})^2)+\omega_{n}^{k+\frac{n}{2}}FR((\omega_{n}^{k+\frac{n}{2}})^2)\)

\(\displaystyle = FL(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}FR(\omega_{n}^{2k+n})\)

\(\displaystyle = FL(\omega_{n}^{2k})+\omega_{n}^{k+\frac{n}{2}}FR(\omega_{n}^{2k})\)

\(\displaystyle = FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})\)

由此，只要知道 \(FL(x)\) 与 \(FR(x)\) 在 \(\displaystyle \omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) 处的取值，就可以得到 \(F(x)\) 在 \(\displaystyle \omega_{n}^{\frac{n}{2}}, \omega_{n}^{\frac{n}{2}+1}, \cdots, \omega_{n}^{n-1}\) 处的取值

拼起来，我们已经可以得到 \(F(x)\) 在 \(\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}\) 处的所有值了

总结下两个关键式子：

\(\bm{\displaystyle F(\omega_{n}^{k}) = FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})}\)
\(\bm{\displaystyle F(\omega_{n}^{k+\frac{n}{2}})= FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}FR(\omega_{\frac{n}{2}}^{k})}\)

据此递归求解即可，时间复杂度 \(O(n \log n)\)

FFT 实现细节

求单位根

小技巧：使用 acos(-1) 可以取到精确的 \(\pi\) 值

由三角函数相关知识，\(\displaystyle \bm{\omega_{n}^{1} = (\cos(\frac{2\pi}{n}), \sin(\frac{2\pi}{n}))}\)，从 \((1, 0)\) 开始一直乘即可

IDFT 的实现

结论：将 DFT 过程中的 \(\omega_{n}^{i}\) 换为 \(\omega_{n}^{-i}\)，做完后除以 \(n\) 即得到系数表达

证明：设点值序列为 \(G(x)\)，则有 \(\displaystyle G(k) = \sum_{i=0}^{n-1} F_i(\omega_{n}^{k})^i\)

现在要证 \(\displaystyle n \times F_k = \sum_{i=0}^{n-1} G_i(\omega_{n}^{-k})^i\)；考虑暴力带入：
\(\displaystyle \sum_{i=0}^{n-1} G_i(\omega_{n}^{-k})^i = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} F_j(\omega_{n}^{i})^j (\omega_{n}^{-k})^i = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} F_j(\omega_{n}^{ij-ik}) = \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} F_j(\omega_{n}^{i(j-k)})\)

考虑对 \(j\) 与 \(k\) 的大小分类讨论：

\(j = k\)，原式为 \(\displaystyle \sum_{i=0}^{n-1} F_k(\omega_{n}^{0}) = n \times F_k\)

\(j \ne k\)，对每个 \(j\)，不妨设 \(t = j-k\)，原式为 \(\displaystyle \sum_{i=0}^{n-1} F_{t+k}(\omega_{n}^{it})\)

由等比数列求和，\(\displaystyle \sum_{i=0}^{n-1} \omega_{n}^{it} = \frac{\omega_{n}^{nt}-\omega_{n}^{0}}{\omega_{n}^{t}-1} = \frac{\omega_{n}^{0}-\omega_{n}^{0}}{\omega_{n}^{t}-1} = 0\)

那么每个 \(j \ne k\) 的 \(j\) 都没有贡献

综上，总贡献为 \(n \times F_k\)，证毕

void fft(cp *f, int l, bool flg){ // 0->IDFT, 1->DFT
	if (l == 1) return;
	cp *fl = f, *fr = f+l/2;
	for (int i=0; i<l; i++) tmp[i] = f[i];
	for (int i=0; i<l/2; i++) fl[i] = tmp[(i<<1)], fr[i] = tmp[(i<<1|1)];
	fft(fl, l/2, flg); fft(fr, l/2, flg); cp w1(cos(2.0*Pi/l), sin(2.0*Pi/l)), w(1, 0);
	if (flg == 0) w1.y = -w1.y;
	for (int i=0; i<l/2; i++){ tmp[i] = fl[i]+w*fr[i]; tmp[i+l/2] = fl[i]-w*fr[i]; w = w*w1; }
	for (int i=0; i<l; i++) f[i] = tmp[i];
}

蝶形运算优化

每次分治，我们将多项式 \(F(x)\) 按照次数奇偶分成左右两半，使用了大量内存移动

若我们一开始就知道 \(F_i\) 最终被分到哪个地方，就可以省去分开的过程，直接倍增合并，达到优化常数的目的

不妨设 \(p(n, k)\) 表示项数为 \(2^n\) 时，\(k\) 次项最终被分到哪里：

若 \(n = 0\)，显然 \(p(n, k) = k\)，原地不动
若 \(n > 0\)，有 \(p(n, k) = ((k\ \texttt{&}\ 1)\ \texttt{<<}\ (n-1))+p(n-1, k\ \texttt{>>}\ 1)\)，奇数项放右边，偶数项放左边

展开可得：\(\displaystyle p(n, k) = \sum_{i=0}^{n-1} (((k\ \texttt{>>}\ i)\ \texttt{&}\ 1)\ \texttt{<<}\ (n-i+1))\)

这显然是将 \(k\) 二进制反转

于是我们预处理数组 \(tr_i\) 表示 \(i\) 的二进制反转，在 FFT 前对每组 \((i, tr_i)\) 对应的系数交换一次即可

下面考虑如何预处理 \(tr_i\)；我们将数 \(i\) 看作最低位与其他位两部分组成

\((tr_{(i\ \texttt{>>}\ 1)})\ \texttt{>>}\ 1\) 可求出其他位的二进制反转；多右移一次是因为反转前 \(i\ \texttt{>>}\ 1\) 时在最高位多带进来了个 \(0\)
若 \(i\ \texttt{&}\ 1\) 为 \(1\)，反转后在最高位补上 \(1\) 即可

预处理代码为 for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0));

合并过程优化

蝶形运算优化后，注意到项数为 \(2^n\) 时，对于 \(k\) 次项：

\(\displaystyle FL(\omega_{\frac{n}{2}}^{k})\) 存储在数组下标为 \(k\) 的位置，\(\displaystyle F(\omega_{n}^{k})\) 应该存储在数组下标为 \(k\) 的位置
\(\displaystyle FR(\omega_{\frac{n}{2}}^{k})\) 存储在数组下标为 \(k+\frac{n}{2}\) 的位置，\(\displaystyle F(\omega_{n}^{k})\) 应该存储在数组下标为 \(k+\frac{n}{2}\) 的位置

因此新开临时数组去记录是非常浪费的，只需要一个临时变量即可

关键代码 (可以结合后面的完整代码理解)：

cp tmp = w*f[k+nl];
f[k+nl] = f[k]-tmp;
f[k] = f[k]+tmp; w = w*w1;

三次变两次

注意到复数有 \(x\) 与 \(y\) 两维，可以同时存储两份信息；我们考虑利用这一性质

记复多项式 \(H(x) = F(x)+G(x)i\)，则有：

\(\displaystyle H^2(x) = (F(x)+G(x)i)^2 = F^2(x)-G^2(x)+2F(x)G(x)i\)

那么只需要一次 DFT，一次 IDFT 即可，最终记得多除 \(2\)

最终代码

P3803 【模板】多项式乘法（FFT）

#include<bits/stdc++.h>
using namespace std;
#define db double
const db Pi = acos(-1);
const int MAXN = 4e6+5;

int n, m, l;
int tr[MAXN];

struct cp{
	db x, y; cp(db xx = 0, db yy = 0){ x = xx; y = yy; }
	cp operator + (const cp& b){ return cp(x+b.x, y+b.y); }
	cp operator - (const cp& b){ return cp(x-b.x, y-b.y); }
	cp operator * (const cp& b){ return cp(x*b.x-y*b.y, x*b.y+y*b.x); }
	cp operator / (const cp& b){ db t = (b.x*b.x+b.y*b.y); return cp((x*b.x+y*b.y)/t, (y*b.x-x*b.y)/t); }
}f[MAXN];

void fft(cp *f, int l, bool flg){ // 0->IDFT, 1->DFT
	for (int i=0; i<l; i++) if (i < tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){
		int nl = (i>>1);
		cp w1(cos(2*Pi/i), sin(2*Pi/i));
		if (flg == 0) w1.y = -w1.y;
		for (int j=0; j<l; j+=i){
			cp w(1, 0);
			for (int k=j; k<j+nl; k++){
				cp tmp = w*f[k+nl];
				f[k+nl] = f[k]-tmp;
				f[k] = f[k]+tmp; w = w*w1;
			}
		}
	}
}

int main(){
	scanf("%d%d", &n, &m);
	for (int i=0; i<=n; i++) scanf("%lf", &f[i].x);
	for (int i=0; i<=m; i++) scanf("%lf", &f[i].y);
	for (l=1; l<=n+m; l<<=1); 
	for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0));
	fft(f, l, 1); for (int i=0; i<l; i++) f[i] = f[i]*f[i]; fft(f, l, 0);
	for (int i=0; i<=n+m; i++) printf("%d ", (int)(f[i].y/l/2+0.49));
	return 0;
}

2. NTT

NTT 原理

FFT 需要依赖复数单位根，容易引发精度误差

大多数问题都在模意义下进行，因此，我们希望为 FFT 找一个模意义下的替代；于是考虑引入原根 \(g\) 代替单位根

不过，对于模数 \(p\) 的原根 \(g\)，其阶数为 \(\varphi(p)=p-1\)；现在对于 \(n = 2^k\ (k \in \mathbb{N^*})\)，我们希望找到阶数为 \(n\) 的数 \(x\)，以此保证 \(x^0, x^1, \cdots, x^{n-1} (\bmod p)\) 互不相同；显然直接用 \(g\) 是不行的

注意到 \(g^k\) 的阶数为 \(\displaystyle \frac{p-1}{\gcd(k, p-1)}\)，取 \(\displaystyle k = \frac{p-1}{n}\) 即可

对 \(g^k\) 阶数的证明：可将 \(g^1, g^2, \cdots, g^{p-1}\) 看成环，每次乘 \(g^k\) 相当于环上走 \(k\) 步

\(k\ \nmid\ (p-1)\)，显然可以走到环上所有位置，阶数为 \(\displaystyle p-1 = \frac{p-1}{\gcd(k, p-1)}\)

\(k\ |\ (p-1)\)，不妨设 \(\displaystyle k = \frac{p-1}{i}\)，则共能走到 \(0, k, 2k, 3k, \cdots, (i-1)k\ (\bmod p)\) 共 \(i\) 个位置，阶数为 \(\displaystyle i = \frac{p-1}{k} = \frac{p-1}{\gcd(k, p-1)}\)

记 \(\displaystyle t = g^{\frac{p-1}{n}}\) (或 \(t_n\))，下面说明 \(t\) 具有用到的所有单位根的性质：

\(\omega_{n}^{k} = \omega_{n}^{k \bmod n}\)：在小环上转圈与在单位圆上转圈类似，多转 \(n\) 次就转回来了，易得 \(t^{k} = t^{k \bmod n}\)；以下省略指数上取模
\(\omega_{n}^{k} = (\omega_{n}^{1})^k\)，更进一步，\(\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i+j}\)：显然 \(t^i \times t^j = t^{i+j}\)
\(\omega_{n}^{0}, \omega_{n}^{1}, \cdots, \omega_{n}^{n-1}\) 互不相同：小环大小为 \(n\)，显然 \(t^{0}, t^{1}, \cdots, t^{n-1}\) 互不相同
\(\displaystyle \omega_{n}^{k+\frac{n}{2}} = -\omega_{n}^{k}\)：
- 即证 \(\displaystyle t^{k+\frac{n}{2}} \equiv -t^{k}\ (\bmod p)\)，移项，即 \(\displaystyle t^{k}(t^{\frac{n}{2}}+1) \equiv 0\ (\bmod p)\)；下面证 \(\displaystyle t^{\frac{n}{2}} \equiv -1\ (\bmod p)\)
- 显然有 \(t^n \equiv 1\ (\bmod p)\) (带入即证)，即 \(\displaystyle (t^{\frac{n}{2}})^2 \equiv 1\ (\bmod p)\)；那么 \(\displaystyle t^{\frac{n}{2}} \bmod p\) 的取值只能为 \(\pm 1\)
- 若 \(\displaystyle t^{\frac{n}{2}} \equiv 1\ (\bmod p)\)，\(t\) 的阶数就不是 \(n\)，而是 \(\displaystyle \frac{n}{2}\) 了，显然矛盾；那么 \(\displaystyle t^{\frac{n}{2}} \equiv -1\ (\bmod p)\)，证毕
\(\omega_{2n}^{2k} = \omega_{n}^{k}\)：直接带入，\(\displaystyle t_{2n} = (g^{\frac{p-1}{2n}})^2 = g^{\frac{p-1}{n}} = t_{n}\)

NTT 实现细节

我们直接将 FFT 中的单位根 \(\omega_{n}^{k}\) 换成 \(t^k\) 即可，求 \(t^1\) 可以直接快速幂，IDFT 中遇到的 \(t^{-1}\) 即为逆元

NTT 要求 \(n\ |\ (p-1)\)，换句话说，\(p-1\) 中需要含有大量 \(2\) 的幂；同时可能需要记住常见模数的原根

P3803 NTT 代码：

#include<bits/stdc++.h>
using namespace std;
#define db double
#define ll long long
const db Pi = acos(-1);
const int MAXN = 4e6+5;
const ll MOD = 998244353, G = 3, INVG = 332748118;

int n, m, l;
int tr[MAXN], f[MAXN], g[MAXN];

ll QuickPow(ll a, ll b){
	ll res = 1;
	while (b){
		if (b&1) res = (res*a)%MOD;
		a = (a*a)%MOD;
		b >>= 1;
	} return res;
}

void ntt(int *f, int l, bool flg){ // 0->IDFT, 1->DFT
	for (int i=0; i<l; i++) if (i < tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){
		int nl = (i>>1), w1 = QuickPow((flg == 1 ? G : INVG), (MOD-1)/i);
		for (int j=0; j<l; j+=i){
			ll w = 1;
			for (int k=j; k<j+nl; k++){
				int tmp = (w*f[k+nl])%MOD;
				f[k+nl] = f[k]-tmp; if (f[k+nl] < 0) f[k+nl] += MOD;
				f[k] = f[k]+tmp; if (f[k] >= MOD) f[k] -= MOD; 
				w = (w*w1)%MOD;
			}
		}
	}
}

int main(){
	scanf("%d%d", &n, &m);
	for (int i=0; i<=n; i++) scanf("%d", &f[i]);
	for (int i=0; i<=m; i++) scanf("%d", &g[i]);
	for (l=1; l<=n+m; l<<=1); int inv = QuickPow(l, MOD-2);
	for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0));
	ntt(f, l, 1); ntt(g, l, 1); for (int i=0; i<l; i++) f[i] = (1ll*f[i]*g[i])%MOD; ntt(f, l, 0);
	for (int i=0; i<=n+m; i++) printf("%d ", (1ll*f[i]*inv)%MOD);
	return 0;
}

这里也放个封装版本：

void pre(int l){ for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0)); }
void ntt(int *g, int l, bool flg){ // 0->IDFT, 1->DFT, l>n
	pre(l); static int f[MAXN] = {0};
	for (int i=0; i<l; i++) f[i] = g[i]; for (int i=0; i<l; i++) if (i<tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){ int nl = (i>>1), w1 = QuickPow((flg?G:INVG), (MOD-1)/i);
		for (int j=0; j<l; j+=i){ ll w = 1;
			for (int k=j; k<j+nl; k++){
				int tmp = (w*f[k+nl])%MOD;
				f[k+nl] = f[k]-tmp; if (f[k+nl]<0) f[k+nl]+=MOD;
				f[k] = f[k]+tmp; if (f[k]>=MOD) f[k]-=MOD; w = (w*w1)%MOD;}}}
	if (flg == 0){ int inv = QuickPow(l, MOD-2); for (int i=0; i<l; i++) g[i] = (1ll*f[i]*inv)%MOD; }
	else{ for (int i=0; i<l; i++) g[i] = f[i]; }
	for (int i=0; i<(l<<1); i++) f[i] = 0;           
}
void mulp(int *f, int *g, int l){ for (int i=0; i<l; i++) f[i] = (1ll*f[i]*g[i])%MOD; } // l > n
void times(int *f, int *g, int n, int m){ // max{len} = n, mod x^m
	static int tmp[MAXN] = {0};
	int l = 1; for (l=1; l<=n+n; l<<=1); for (int i=0; i<l; i++) tmp[i] = g[i];
	ntt(f, l, 1); ntt(tmp, l, 1); mulp(f, tmp, l); ntt(f, l, 0);
	for (int i=0; i<l; i++) tmp[i] = 0; for (int i=m; i<l; i++) f[i] = 0;
}

函数多了很容易搞混，包括后面，实现时千万注意传的参是次数还是 \(n = 2^k\ (k \in \mathbb{N^*})\)！

3. 多项式全家桶

多项式四则运算

多项式加减法

多项式的加减定义为：

\(\displaystyle F(x)+G(x) = \sum_{i=0} (F_i+G_i)x^i\)
\(\displaystyle F(x)-G(x) = \sum_{i=0} (F_i-G_i)x^i\)

系数表达直接相加 / 相减即可，可以 \(O(n)\) 算出

多项式乘法

使用上文所述 FFT 或 NTT 即可，可以 \(O(n \log n)\) 算出

多项式乘法逆元

当 \(F(x)G(x) = 1\) 时，我们称 \(F(x)\) 与 \(G(x)\) 互为乘法逆元；其可帮助我们将(整)除转化为乘

考虑次数倍增，不断扩大次数上界，直到求出完整的逆元

不妨设 \(R_{\triangle}(x)\) 为 \(F(x)\) 的前 \(\displaystyle \frac{n}{2}\) 位逆元，考虑如何求出 \(R(x)\) 表示 \(F(x)\) 的前 \(n\) 位逆元

显然有：\(\displaystyle R(x)-R_{\triangle}(x) \equiv 0\ (\bmod x^{\frac{n}{2}})\)；考虑平方求得 \(\bmod x^n\) 下的逆元：

\(\displaystyle (R(x)-R_{\triangle}(x))^2 \equiv 0\ (\bmod x^{n})\)

\(\displaystyle R^2(x)-2R(x)R_{\triangle}(x)+R_{\triangle}^2(x) \equiv 0\ (\bmod x^{n})\)，乘个 \(F(x)\) 消去 \(R(x)\)：

\(\displaystyle F(x)(R^2(x)-2R(x)R_{\triangle}(x)+R_{\triangle}^2(x)) \equiv 0\ (\bmod x^{n})\)

\(\displaystyle R(x)-2R_{\triangle}(x)+R_{\triangle}^2(x)F(x) \equiv 0\ (\bmod x^{n})\)

\(\displaystyle R(x) \equiv 2R_{\triangle}(x)-R_{\triangle}^2(x)F(x)\ (\bmod x^{n})\)

根据上式倍增即可；边界条件为 \(\bmod x^1\) 时就是求数字逆元

void invp(int *f, int n){ // len = n, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}, tmp[MAXN] = {0}; w[0] = QuickPow(f[0], MOD-2);
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; }
		for (int j=0; j<i; j++) tmp[j] = f[j];
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(tmp, (i<<1), 1); mulp(w, tmp, (i<<1)); ntt(w, (i<<1), 0); 
        // 三个多项式相乘，次数超过 i，开到 (i<<1) 避免循环卷积溢出破坏前面几位
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = r[j]-w[j]; if (w[j] < 0) w[j] += MOD; }
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = tmp[i] = 0;
}

多项式带余除法

给定 \(n\) 次多项式 \(F(x)\) 与 \(m\) 次多项式 \(G(x)\)，令 \(F(x) = Q(x)G(x)+R(x)\)，求 \(n-m\) 次商多项式 \(Q(x)\) 与 \(m-1\) 次余多项式 \(R(x)\)

我们希望利用 \(\bmod x^k\) 的方式消去 \(R(x)\)；但这种方法只适用于消高次项，不适用于消低次项，于是我们考虑将多项式反转

定义 \(n\) 次多项式 \(F(x)\) 的反转为 \(\displaystyle F^{T}(x) = x^nF(\frac{1}{x})\)，容易发现即为将其系数前后反转

例：\(F(x) = 3+x+4x^2\)，则有 \(\displaystyle F^T(x) = x^2F(\frac{1}{x}) = 3x^2+x^2 \cdot \frac{1}{x}+x^2 \cdot \frac{4}{x^2} = 4+x+3x^2\)

那么由 \(F(x) = Q(x)G(x)+R(x)\)，考虑同乘 \(x^n\) 转化为反转形式：

\(\displaystyle x^nF(\frac{1}{x}) = (x^{n-m}Q(\frac{1}{x}))(x^{m}G(\frac{1}{x}))+x^{n-m+1}(x^{m-1}R(\frac{1}{x}))\)

\(\displaystyle F^T(x) = Q^T(x)G^T(x)+x^{n-m+1}R^T(x)\)，此时对 \(x^{n-m+1}\) 取模消去 \(R^T(x)\)：

\(\displaystyle F^T(x) = Q^T(x)G^T(x) (\bmod x^{n-m+1})\)

\(\displaystyle Q^T(x) = F^T(x)(G^T(x))^{-1} (\bmod x^{n-m+1})\)

先全部反转，求逆后相乘可得 \(Q^T(x)\)，反转回来即得 \(Q(x)\)；最后将 \(Q(x)\) 与 \(G(x)\) 乘起来，与 \(F(x)\) 作差即得 \(R(x)\)

反转可以直接用自带的 reverse 函数

void modp(int *f, int *g, int n, int m){ // len_f = n, len_g = m, len_q = n-m, len_r = m-1
	static int g_[MAXN] = {0}, f_[MAXN] = {0}; int l = n-m+1;
	reverse(g, g+m+1); for (int i=0; i<l; i++) g_[i] = g[i]; reverse(g, g+m+1);
	reverse(f, f+n+1); for (int i=0; i<l; i++) f_[i] = f[i]; reverse(f, f+n+1);
	invp(g_, l); times(g_, f_, l, l); reverse(g_, g_+l); times(g, g_, n, n); 
	for (int i=0; i<n; i++){ g[i] = f[i]-g[i]; if (g[i] < 0) g[i] += MOD; }
	for (int i=0; i<l; i++) f[i] = g_[i]; for (int i=m; i<n; i++) g[i] = 0; 
	for (int i=0; i<(l<<1); i++) g_[i] = f_[i] = 0;
}

牛顿迭代

牛顿迭代是一种常用的分析方法，解决形如 "\(G(F(x)) = 0\)，求 \(F(x)\ (\bmod x^n)\)" 的问题

这里不加证明地给出结论：\(\displaystyle \bm{F(x) = (F_{\triangle}(x)-\frac{G(F_{\triangle}(x))}{G'(F_{\triangle}(x))})\ (\bmod x^n)}\) (特别重要)

其中，\(F_{\triangle}(x)\) 表示倍增上一层的求解结果，\(G'(x)\) 表示多项式 \(G(x)\) 求导

\(n\) 次多项式 \(\displaystyle F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\) 求导定义为：\(\displaystyle F'(x) = a_1+2a_2x+3a_3x^2+\cdots+(n-1)a_{n-1}x^{n-1}\)
\(n\) 次多项式 \(\displaystyle F(x) = a_0+a_1x+a_2x^2+\cdots+a_nx^n\) 积分定义为：\(\displaystyle \int F(x)\text{dx} = C+a_0x+\frac{1}{2}a_1x^2+\cdots+\frac{1}{n+1}a_nx^{n+1}\)
实现时需要线性求逆元：
void Invinit(int n){ inv[1] = 1; for (int i=2; i<=n; i++) inv[i] = (1ll*inv[MOD%i]*(MOD-MOD/i))%MOD; }
void derip(int *f, int n){ for (int i=1; i<=n; i++) f[i-1] = (1ll*f[i]*i)%MOD; f[n] = 0; } // len = n, F'(x)
void sump(int *f, int n){ Invinit(n+1); for (int i=n+1; i>=1; i--) f[i] = (1ll*f[i-1]*inv[i])%MOD; f[0] = 0; } // len = n

多项式开根

已知 \(B(x)\)，希望求使 \(A^2(x)-B(x) = 0\) 的 \(A(x)\)

法一：构造 \(G(A(x)) = A^2(x)-B(x)\)，则其导数 \(G'(A(x)) = 2A(x)\) (将 \(B(x)\) 看作常量)，考虑直接套牛顿迭代：

\(\displaystyle A(x) = A_{\triangle}(x)-\frac{G(A_{\triangle}(x))}{G'(A_{\triangle}(x))} = A_{\triangle}(x)-\frac{A_{\triangle}^2(x)-B(x)}{2A_{\triangle}(x)} = \frac{A_{\triangle}^2(x)+B(x)}{2A_{\triangle}(x)}\)

那么直接倍增，求逆+乘法即可

边界条件为：

若 \(B_0 = 1\)，则 \(A_0 = 1\)
反之，需要求解同余方程 \(A_0^2 \equiv B_0 (\bmod p)\)

法二：转化为后文要讲的多项式快速幂，指数视为 \(2\) 关于 \(p\) 的逆元即可

void sqrtp(int *f, int n){ // rem->f, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; } invp(r, i);
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(w, (i<<1), 0);
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = f[j]+w[j]; if (w[j] >= MOD) w[j] -= MOD; } times(w, r, i, i); r[i] = 0;
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}

多项式 ln 与 exp

两者均由麦克劳林级数定义：

\(\displaystyle \ln(F(x)) = -\sum_{i=1} \frac{(1-F(x))^i}{i}\)
\(\displaystyle \exp(F(x)) = \sum_{i=0} \frac{F^i(x)}{i!}\)

性质：\(\ln(F(x))\) 与 \(\exp(F(x))\) 互为逆运算 (重要)

对于 \(\ln(F(x))\)，我们肯定是不用定义式来求的；注意到 \(\displaystyle \ln'(x) = \frac{1}{x}\)，于是对 \(\displaystyle \ln(A(x)) = B(x)\) 同时求导：

\(\displaystyle \ln'(A(x)) = \frac{A'(x)}{A(x)} = B'(x)\) (注意复合函数求导法则)

那么对 \(\displaystyle \frac{A'(x)}{A(x)}\) 再积分一次即得 \(B(x)\)

void lnp(int *f, int n){ // len = n
	static int g[MAXN] = {0};
	for (int i=0; i<=n; i++) g[i] = f[i]; invp(g, n);
	derip(f, n); times(f, g, n, n); sump(f, n-1);
	for (int i=0; i<=n; i++) g[i] = 0;
}

对于 \(\exp(F(x))\)，由 \(\displaystyle e^{A(x)} = e^{\ln(B(x))} = B(x)\)，考虑构造 \(G(B(x)) = \ln(B(x))-A(x)\)，直接牛顿迭代

\(\displaystyle G(B(x)) = B_{\triangle}(x)-\frac{G(B_{\triangle}(x))}{G'(B_{\triangle}(x))} = B_{\triangle}(x)-\frac{\ln(B_{\triangle}(x))-A(x)}{\frac{1}{B_{\triangle}(x)}} = B_{\triangle}(x)(1-\ln(B_{\triangle}(x))+A(x))\)

那么直接倍增，求 \(\ln\) + 乘法即可；边界条件为 \(A_0 = 1\) 时，\(B_0 = 1\)

void exp(int *f, int n){ // len = n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = w[j]; } lnp(r, i);
		for (int j=0; j<i; j++){ r[j] = f[j]-r[j]; if (r[j] < 0) r[j] += MOD; } 
		r[0] = (r[0]+1)%MOD; times(w, r, i, i); r[i] = 0;
 	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}

多项式快速幂

给定 \(A(x)\) 与 \(k\)，求 \(A^k(x)\)；考虑对 \(A^k(x)\) 取 \(\ln\)，将乘法变为加法，最后 \(\exp\) 回去即可

式子为 \(\displaystyle A^k(x) = \exp(k \ln(A(x)))\)；时间复杂度 \(O(n \log n)\)

细节：如果指数 \(k > p\)，应直接对 \(p\) 取模而非 \(p-1\)

证明：将 \(\exp(kF(x))\) 按定义式展开，为 \(\displaystyle \sum_{i=0} \frac{k^iF^i(x)}{i!}\)，显然 \(k^i\) 应对 \(p\) 取模

void powp(int *f, int n, int k){ lnp(f, n); for (int i=0; i<=n; i++) f[i] = (1ll*f[i]*k%MOD+MOD)%MOD; exp(f, n); } // len = n

任意模数 NTT (MTT)

当模数并非 \(r \times 2^k+1\) 的形式时，我们需要精度足够高的算法解决

直接 FFT 或 NTT 的值域太大，我们无法接受

考虑拆系数，将每个系数拆为 \(a_1 \times 2^{15}+a_2\) 的形式，此时值域大概为 \(n \times 2^{15} \times 2^{15} \approx 10^{15}\)

那么原来两个系数 \(c_1 = a_1 \times 2^{15}+a_2\) 与 \(c_2 = b_1 \times 2^{15}+b_2\) 的乘积即为：

\((a_1 \times 2^{15} + a_2)(b_1 \times 2^{15} + b_2) = a_1b_12^{30}+(a_1b_2+a_2b_1)2^{15}+a_2b_2\)

直接进行四次多项式乘法显然会常数爆炸，不过我们可以考虑利用复数的性质同时计算多份信息 (类似三次变两次优化)

设复多项式 \(P = A_1+iA_2\)，\(Q = B_1+iB_2\)，则有 \(PQ = A_1B_1-A_2B_2+(A_1B_2+A_2B_1)i\)

类似共轭，取 \(P = A_1-iA_2\)，则有 \(P'Q = A_1B_1+A_2B_2+(A_1B_2-A_2B_1)i\)

将 \(PQ\) 与 \(P'Q\) 相加，可得 \(A_1B_1\) 与 \(A_1B_2\)

用 \(A_1B_2+A_2B_1\) 减去 \(A_1B_2\) 可得 \(A_2B_1\)，用 \(A_1B_1+A_2B_2\) 减去 \(A_1B_1\) 可得 \(A_2B_2\)

那么总共只需要三次 DFT，两次 IDFT 即可

P4245 【模板】任意模数多项式乘法

#include<bits/stdc++.h>
using namespace std;
#define db long double
#define ll long long
const db Pi = acos(-1);
const int MAXN = 4e6+5;

int n, m, p;
int tr[MAXN], f[MAXN], g[MAXN];

struct cp{
	db x, y; cp(db xx = 0, db yy = 0){ x = xx; y = yy; }
	cp operator + (const cp& b){ return cp(x+b.x, y+b.y); }
	cp operator - (const cp& b){ return cp(x-b.x, y-b.y); }
	cp operator * (const cp& b){ return cp(x*b.x-y*b.y, x*b.y+y*b.x); }
	cp operator / (const cp& b){ db t = (b.x*b.x+b.y*b.y); return cp((x*b.x+y*b.y)/t, (y*b.x-x*b.y)/t); }
};

void pre(int l){ for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0)); }
void fft(cp *f, int l, bool flg){ // 0->IDFT, 1->DFT
	pre(l); for (int i=0; i<l; i++) if (i<tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){ int nl = (i>>1); cp w1(cos(2*Pi/i), sin(2*Pi/i)); if (flg == 0) w1.y = -w1.y;
		for (int j=0; j<l; j+=i){ cp w(1, 0);
			for (int k=j; k<j+nl; k++){
				cp tmp = w*f[k+nl]; f[k+nl] = f[k]-tmp; f[k] = f[k]+tmp; w = w*w1;}}}
	if (flg == 0) for (int i=0; i<l; i++) f[i].x /= l, f[i].y /= l;
}
void mtt(int *f, int *g, int n, int m){
	static cp p1[MAXN], p2[MAXN], q[MAXN]; 
	int l = 0; for (l=1; l<=n+n; l<<=1);
	for (int i=0; i<l; i++){
		p1[i] = cp(f[i]>>15, f[i]&32767);
		p2[i] = cp(f[i]>>15, -(f[i]&32767));
		q[i] = cp(g[i]>>15, g[i]&32767);
	} fft(p1, l, 1); fft(p2, l, 1); fft(q, l, 1);
	for (int i=0; i<l; i++){ p1[i] = p1[i]*q[i]; p2[i] = p2[i]*q[i]; } fft(p1, l, 0); fft(p2, l, 0);
	for (int i=0; i<m; i++){
		ll a = (ll)floor((p1[i].x+p2[i].x)/2+0.49)%p;
		ll b = (ll)floor((p1[i].y+p2[i].y)/2+0.49)%p;
		ll c = (ll)floor((p1[i].y-b)+0.49)%p;
		ll d = (ll)floor((p2[i].x-a)+0.49)%p;
		f[i] = ((((a<<15)+(b+c))<<15)+d)%p;
		if (f[i] < 0) f[i] += p;
	}
}

int main(){
	scanf("%d%d%d", &n, &m, &p);
	for (int i=0; i<=n; i++) scanf("%d", &f[i]);
	for (int i=0; i<=m; i++) scanf("%d", &g[i]); mtt(f, g, max(n, m), n+m+1);
	for (int i=0; i<=n+m; i++) printf("%d ", f[i]);
	return 0;
}

4. 代码打包

没加上 MTT

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MAXN = 4e6+5;
const ll MOD = 998244353, G = 3, INVG = 332748118, INV2 = 499122177;

int n, k;
int tr[MAXN], f[MAXN], g[MAXN], inv[MAXN];

int read(){
	bool f=1;ll x=0;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-') f=!f;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);x=x%MOD;ch=getchar();}
	x=(f?x:-x);return (int)x;
}

ll QuickPow(ll a, ll b){
	ll res = 1;
	while (b){
		if (b&1) res = (res*a)%MOD;
		a = (a*a)%MOD;
		b >>= 1;
	} return res;
}

void Invinit(int n){ inv[1] = 1; for (int i=2; i<=n; i++) inv[i] = (1ll*inv[MOD%i]*(MOD-MOD/i))%MOD; }
void pre(int l){ for (int i=0; i<l; i++) tr[i] = ((tr[(i>>1)]>>1)|((i&1)?(l>>1):0)); }
void ntt(int *g, int l, bool flg){ // 0->IDFT, 1->DFT, l>n
	pre(l); static int f[MAXN] = {0};
	for (int i=0; i<l; i++) f[i] = g[i]; for (int i=0; i<l; i++) if (i<tr[i]) swap(f[i], f[tr[i]]);
	for (int i=2; i<=l; i<<=1){ int nl = (i>>1), w1 = QuickPow((flg?G:INVG), (MOD-1)/i);
		for (int j=0; j<l; j+=i){ ll w = 1;
			for (int k=j; k<j+nl; k++){
				int tmp = (w*f[k+nl])%MOD;
				f[k+nl] = f[k]-tmp; if (f[k+nl]<0) f[k+nl]+=MOD;
				f[k] = f[k]+tmp; if (f[k]>=MOD) f[k]-=MOD; w = (w*w1)%MOD;}}}
	if (flg == 0){ int inv = QuickPow(l, MOD-2); for (int i=0; i<l; i++) g[i] = (1ll*f[i]*inv)%MOD; }
	else{ for (int i=0; i<l; i++) g[i] = f[i]; }
	for (int i=0; i<(l<<1); i++) f[i] = 0;           
}
void mulp(int *f, int *g, int l){ for (int i=0; i<l; i++) f[i] = (1ll*f[i]*g[i])%MOD; } // l > n
void times(int *f, int *g, int n, int m){ // max{len} = n, mod x^m
	static int tmp[MAXN] = {0};
	int l = 1; for (l=1; l<=n+n; l<<=1); for (int i=0; i<l; i++) tmp[i] = g[i];
	ntt(f, l, 1); ntt(tmp, l, 1); mulp(f, tmp, l); ntt(f, l, 0);
	for (int i=0; i<l; i++) tmp[i] = 0; for (int i=m; i<l; i++) f[i] = 0;
}
void invp(int *f, int n){ // len = n, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}, tmp[MAXN] = {0}; w[0] = QuickPow(f[0], MOD-2);
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; }
		for (int j=0; j<i; j++) tmp[j] = f[j];
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(tmp, (i<<1), 1); mulp(w, tmp, (i<<1)); ntt(w, (i<<1), 0); 
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = r[j]-w[j]; if (w[j] < 0) w[j] += MOD; }
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = tmp[i] = 0;
}
void sqrtp(int *f, int n){ // rem->f, mod x^n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = (w[j]<<1); if (r[j] >= MOD) r[j] -= MOD; } invp(r, i);
		ntt(w, (i<<1), 1); mulp(w, w, (i<<1)); ntt(w, (i<<1), 0);
		for (int j=i; j<(i<<1); j++) w[j] = 0;
		for (int j=0; j<i; j++){ w[j] = f[j]+w[j]; if (w[j] >= MOD) w[j] -= MOD; } times(w, r, i, i); r[i] = 0;
	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}
void modp(int *f, int *g, int n, int m){ // len_f = n, len_g = m, len_q = n-m, len_r = m-1
	static int g_[MAXN] = {0}, f_[MAXN] = {0}; int l = n-m+1;
	reverse(g, g+m+1); for (int i=0; i<l; i++) g_[i] = g[i]; reverse(g, g+m+1);
	reverse(f, f+n+1); for (int i=0; i<l; i++) f_[i] = f[i]; reverse(f, f+n+1);
	invp(g_, l); times(g_, f_, l, l); reverse(g_, g_+l); times(g, g_, n, n); 
	for (int i=0; i<n; i++){ g[i] = f[i]-g[i]; if (g[i] < 0) g[i] += MOD; }
	for (int i=0; i<l; i++) f[i] = g_[i]; for (int i=m; i<n; i++) g[i] = 0; 
	for (int i=0; i<(l<<1); i++) g_[i] = f_[i] = 0;
}
void derip(int *f, int n){ for (int i=1; i<=n; i++) f[i-1] = (1ll*f[i]*i)%MOD; f[n] = 0; } // len = n, F'(x)
void sump(int *f, int n){ Invinit(n+1); for (int i=n+1; i>=1; i--) f[i] = (1ll*f[i-1]*inv[i])%MOD; f[0] = 0; } // len = n
void lnp(int *f, int n){ // len = n
	static int g[MAXN] = {0};
	for (int i=0; i<=n; i++) g[i] = f[i]; invp(g, n);
	derip(f, n); times(f, g, n, n); sump(f, n-1);
	for (int i=0; i<=n; i++) g[i] = 0;
}
void exp(int *f, int n){ // len = n
	int l = 0; for (l=1; l<=n; l<<=1);
	static int r[MAXN] = {0}, w[MAXN] = {0}; w[0] = 1;
	for (int i=2; i<=l; i<<=1){
		for (int j=0; j<(i>>1); j++){ r[j] = w[j]; } lnp(r, i);
		for (int j=0; j<i; j++){ r[j] = f[j]-r[j]; if (r[j] < 0) r[j] += MOD; } 
		r[0] = (r[0]+1)%MOD; times(w, r, i, i); r[i] = 0;
 	} for (int i=0; i<=n; i++) f[i] = w[i];
	for (int i=0; i<(l<<1); i++) r[i] = w[i] = 0;
}
void powp(int *f, int n, int k){ lnp(f, n); for (int i=0; i<=n; i++) f[i] = (1ll*f[i]*k%MOD+MOD)%MOD; exp(f, n); } // len = n

int main(){		
	return 0;
}

posted @ 2025-05-29 15:11 lzlqwq 阅读(32) 评论(0) 收藏举报

刷新页面返回顶部

lzlqwq

多项式基础 - 学习笔记

多项式基础 - 学习笔记

1. FFT

前置知识：复数

基础定义

四则运算

辐角与辐角主值

单位根

DFT 与 IDFT 思想

FFT 原理

FFT 实现细节

求单位根

IDFT 的实现

蝶形运算优化

合并过程优化

三次变两次

最终代码

2. NTT

NTT 原理

NTT 实现细节

3. 多项式全家桶

多项式四则运算

多项式加减法

多项式乘法

多项式乘法逆元

多项式带余除法

牛顿迭代

多项式开根

多项式 ln 与 exp

多项式快速幂

任意模数 NTT (MTT)

4. 代码打包

公告