526互联

多项式（Poly）笔记

发布时间 2023-12-23 09:02:24作者: Lcy++

开头先扔板子：多项式板子们

定义

多项式（polynomial）是形如 \(P(x) = \sum \limits_{i = 0}^{n} a_i x ^ i\) 的代数表达式。其中 \(x\) 是一个不定元。

\(\partial(P(x))\) 称为这个多项式的次数。

多项式的基本运算

  • 多项式的加减法

\[A(x) = \sum_{i = 0}^{n} a_i x ^ i, B(x) = \sum_{i = 0}^{m} b_i x ^ i \]

\[A(x) \pm B(x) = \sum_{i = 0}^{\max(n, m)} (a_i \pm b_i) x ^ i \]

  • 多项式的乘法

\[A(x) \times B(x) = \sum_{k = 0}^{nm} \sum_{i + j = k} a_i b_j x ^ k \]

  • 多项式除法

这里讨论多项式的带余除法。

可以证明，一定存在唯一的 \(q(x), r(x)\) 满足 \(A(x) = q(x) B(x) + r(x)\)，且 \(\partial(r(x)) < \partial(B(x))\)。

\(q(x)\) 称为 \(B(x)\) 除 \(A(x)\) 的商式，\(r(x)\) 称为 \(B(x)\) 除 \(A(x)\) 的余式。记作：

\[A(x) \equiv r(x) (\bmod \ B(x)) \]

特别的，若 \(r(x) = 0\)，则称 \(B(x)\) 整除 \(A(x)\)，\(A(x)\) 是 \(B(x)\) 的一个倍式，\(B(x)\) 是 \(A(x)\) 的一个因式。记作 \(B(x) | A(x)\)。

有关多项式的引理

  • 对于 \(n + 1\) 个点可以唯一确定一个 \(n\) 次多项式。

  • 如果 \(f(x), g(x)\) 都是不超过 \(n\) 次的多项式，且它对于 \(n +1\) 个不同的数 \(\alpha_1, \alpha_2 \cdots \alpha_n\) 有相同的值，即 \(\forall i \in [1, n + 1], i \in \mathbb{Z}, f(\alpha_i) = g(\alpha_i)\)。则 \(f(x) = g(x)\)。

多项式的点值表示

如果选取 \(n + 1\) 个不同的数 \(x_0, x_1, \cdots x_n\) 对多项式进行求值，得到 \(A(x_0), A(x_1) \cdots A(x_n)\)，那么就称 \(\{(x_i, A(x_i)) \ | \ 0 \le i \le n, i \in \mathbb{Z} \}\) 为 \(A(x)\) 的点值表示。

快速傅里叶变换（FFT）

快速傅里叶变换是借助单位根进行求值和插值，从而快速进行多项式乘法的算法。

单位根

将复平面上的单位圆均分成 \(n\) 份，从 \(x\) 轴数，第 \(i\) 条线与单位圆的交点称为 \(i\) 次单位根，记作 \(\omega_{n}^{i}\)。

根据定义，可以得到：\(\omega_{n}^{1} = i \sin \alpha +\cos \alpha, \alpha = \frac{2 \pi}{n}\)。

根据欧拉恒等式，可以得到：\(\omega_{n}^{1} = e ^ {\frac{2 \pi i}{n}}\)。

由此那么可以得到下面的性质：

  • \(\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i + j}\)

  • \(\omega_{n}^{i + \frac{n}{2}} = - \omega_{n}^{i}\)

  • \(\omega_{n}^{i + n} = \omega_{n}^{i}\)

离散傅里叶变换（DFT）

离散傅里叶变换，是将 \(\omega_{n}^{k}, 0 \le k < n\) 代入 \(f(x)\) 和 \(g(x)\) 中求值的过程。

对于朴素的方法，每次代入一个单位根，然后用 \(O(n)\) 的复杂度计算函数值。时间复杂度 \(O(n ^ 2)\)。

离散傅里叶变换利用了单位根的性质巧妙优化了这个过程。离散傅里叶变换过程如下：

首先将 \(f(x)\) 根据次数的奇偶性拆成两部分，分别分为：

\[\begin{cases} e(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i} x ^ {2i} = a_0 + a_2 x^2 + a_4 x^4 \cdots a_{n - 2} x^{n - 2} \\ o(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i + 1} x ^ {2i + 1} = a_1 x + a_3 x^3 + a_5 x^5 \cdots a_{n - 1} x^{n - 1} \end{cases} \]

设

\[\begin{cases} e'(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i} x ^ i = a_0 + a_2 x + a_4 x^2 \cdots a_{n - 2} x^{\frac{n - 2}{2}} \\ o'(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i + 1} x ^ {i} = a_1 + a_3 x^1 + a_5 x^2 \cdots a_{n - 1} x^{\frac{n - 2}{2}} \end{cases} \]

则 \(f(x) = e'(x ^ 2) + x o'(x ^ 2)\)。

将 \(\omega_{n}^{k}\) 代入得到：

\[\begin{cases} f(\omega_{n}^{k}) = e'((\omega_{n}^{k}) ^ 2) + \omega_{n}^{k} o'((\omega_{n}^{k}) ^ 2) = e'(\omega_{n}^{2k}) + \omega_{n}^{k} o'(\omega_{n}^{2k}) \\\\ f(\omega_{n}^{k + \frac{n}{2}}) = e'((\omega_{n}^{k+ \frac{n}{2}}) ^ 2) + \omega_{n}^{k+ \frac{n}{2}} o'((\omega_{n}^{k+ \frac{n}{2}}) ^ 2) = e'(\omega_{n}^{2k}) - \omega_{n}^{k} o'(\omega_{n}^{2k}) \end{cases} \]

然后你发现，\(f(\omega_{n}^{k})\) 和 \(f(\omega_{n}^{k + \frac{n}{2}})\) 仅仅差了一个符号！！！

所以只需要求出 \(e'(x)\) 和 \(o'(x)\) 对 \(\omega^{k}_{n}\)（\(0 \le k \le \frac{n}{2}\)）上的取值，就可以推出 \(f(x)\) 在两倍点数上的取值。

每次问题规模缩小一半，因此时间复杂度 \(O(n \log n)\)。

离散傅里叶逆变换（IDFT）

假设对于两个多项式都得到了 \(t = n + m - 1\) 个点值，设为 \(\{(x_i, A(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}, \{(x_i, B(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}\)。

那么可以知道，多项式 \(C(x) = A(x) \times B(x)\) 的点值表示一定为：

\[\{(x_i, A(x_i) B(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \} \]

现在，只需要将这 \(t\) 个点插值回去，就可以得到 \(A(x)B(x)\) 了。

先设这 \(t\) 个点值分别是：\(\{(x_i, v_i) \ | \ 0 \le i < t, i \in \mathbb{Z} \}\)，设最后的多项式为 \(C(x) = \sum \limits_{i = 0}^{n + m - 2} c_i x^i\)，这里直接给出结论：

\[c_k = \dfrac{1}{n} \sum_{i = 0}^{t - 1} v_i \omega_{t}^{-ki} \]

由此可见，IDFT 和 DFT 仅仅有一个负号的区别。只要将所有的单位根从 \(\omega_{n}^{k}\) 变成 $ - \omega_{n}^{k}$ 即可。

void FFT(cp a[], int n, int op) {
	if (n == 1) return;
	cp a1[n], a2[n];
	rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
	FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
	cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
	cp Wk = {1, 0};
	rop(i, 0, n >> 1) {
		a[i] = a1[i] + Wk * a2[i];
		a[i + (n >> 1)] = a1[i] - Wk * a2[i];
		Wk = Wk * Wn;
	}
}
void FFT(cp a[], cp b[], int n, int m) {
	m = n + m; n = 1;
	while (n <= m) n <<= 1;
	FFT(a, n, 1); FFT(b, n, 1);
	rop(i, 0, n) a[i] = a[i] * b[i];
	FFT(a, n, -1);
	rep(i, 0, m) a[i].x = a[i].x / n;
}

FFT 优化

  • 三次变两次优化

原本的朴素 FFT，将 \(\{a\}, \{b\}\) 两个序列分别求值，乘起来再 IDFT 插值一下，一共跑了三次 FFT。这是不好的。

三次变两次优化是这样的：将原序列合并成一个复数：\(\{a_i + b_i \times i\}\)。做一遍 DFT 把求出的每个函数值平方。因为 \((a + bi) ^ 2 = (a ^ 2 - b ^ 2) + (2abi)\)。因此把虚部取出来以后除以 \(2\) 就是答案。

void FFT(cp a[], int n, int op) {
	if (n == 1) return;
	cp a1[n], a2[n];
	rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
	FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
	cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
	cp Wk = {1, 0};
	rop(i, 0, n >> 1) {
		a[i] = a1[i] + a2[i] * Wk;
		a[i + (n >> 1)] = a1[i] - a2[i] * Wk;
		Wk = Wk * Wn;
	}
}
int main() {
	read(n, m);
	rep(i, 0, n) scanf("%lf", &a[i].x);
	rep(i, 0, m) scanf("%lf", &a[i].y);
	m = n + m; n = 1;
	while (n <= m) n <<= 1;
	FFT(a, n, 1);
	rop(i, 0, n) a[i] = a[i] * a[i];
	FFT(a, n, -1);
	rep(i, 0, m) printf("%d ", (int)(a[i].y / 2 / n + 0.5));
}
  • 蝴蝶变换优化

后面再补吧。其实本人感觉这个优化不是那么必要，因为三次变两次实在太快了。

FFT 例题

  1. A * B problem(plus)

可以设 \(x = 10\)，把 \(a\) 写成 \(A(x) = \sum \limits_{i = 0}^{n} a_i x^i\) 的形式（\(n = \log_{10} a\)）。同理可以把 \(b\) 转化为多项式 \(B(x)\)。

这样求两个数相乘就是求 \(A(x) \times B(x)\) 啊。

所以直接 \(O(n \log n)\) 做完了。

  1. P3338 [ZJOI2014] 力

给出 \(n\) 个数 \(q_1,q_2, \dots q_n\)，定义

\[F_j~=~\sum_{i = 1}^{j - 1} \frac{q_i \times q_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i \times q_j}{(i - j)^2} \]

\[E_i~=~\frac{F_i}{q_i} \]

对 \(1 \leq i \leq n\)，求 \(E_i\) 的值。

首先发现这个除以 \(q_i\) 就是没用。所以可以化简成：

\[E_j = \sum_{i = 1}^{j - 1} \frac{q_i}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i}{(i - j)^2} \]

先看前面这个式子。答案就是：

\[(j - 1) ^ 2 q_1 + (j - 2) ^ 2 q_2 + (j - 3) ^ 2 q_3 \cdots \]

设 \(f(x) = \sum q_i x ^ i, g(x) = \dfrac{1}{i ^ 2} x ^ i\)。可以发现，\(E_j' = (f \times g)[j]\)

再看后面这一块的式子。我们把 \(f(x)\) 的系数翻转，变成 \(f'(x) = \sum q_{n - j + 1} x ^ j\)。那么可以发现 \(E_{j}'' = (f' \times g)[n - j + 1]\)。

跑两次 FFT 就完事了。

  1. P3723 [AH2017/HNOI2017] 礼物

首先发现加减相对于两个手环是对称的。因此可以把对一个手环的减法转化成对另一个手环的加法。这样可以假设全是在第一个手环上执行的加减操作。

第一个手环执行了加 \(c\) 的操作，且旋转过之后的序列为 \([x_1, x_2 \cdots x_n]\)，第二个手环为 \([y_1, y_2 \cdots y_n]\)。计算差异值并化简，可以得到差异值是：

\(\sum x ^ 2 + \sum y ^ 2 + c ^ 2 n + 2c(\sum x - \sum y) - 2 \sum xy\)

可以发现，这个序列只有最后一项是不定的。

因此将 \(y\) 序列翻转后再复制一倍，与 \(x\) 卷积，答案就是卷积后序列的 \(n + 1 \sim 2n\) 项系数的 \(\max\)。

直接暴力枚举 \(c\)，加上前面依托就行了。

\(\text{AC code}\)

快速数论变换（NTT）

快速数论变换就是基于原根的快速傅里叶变换。

首先考虑快速傅里叶变换用到了单位根的什么性质。

  1. \(\omega_{n}^{k}, 0 \le k < n\) 互不相同。

  2. \(\omega_{n}^{k} = \omega_{n}^{k + \frac{n}{2}}\)。

  3. \(\omega_{n}^{k} = \omega_{2n}^{2k}\)。

数论中，原根恰好满足这些性质。

对于一个素数的原根 \(g\)，设 \(g_n = g ^ {\frac{p - 1}{n}}\)。那么：

\[g_n ^ n = g ^ {p - 1} \equiv 1(\bmod \ p) \]

\[g_n ^ {\frac{n}{2}} = g ^ {\frac{p - 1}{2}} \equiv - 1(\bmod ~ p) \]

\[g ^ \alpha + g ^ \beta = g ^ {\alpha + \beta} \]

\[g_{\alpha n}^{\alpha k} = g_{n}^{k} \]

我们发现它满足 \(\omega_{n}^{k}\) 的全部性质！

因此，只需要用 \(g_{n}^k = g_{}^{\frac{p - 1}{n} k}\) 带替全部的 \(\omega_{n}^{k}\) 就可以了。

tips：对于一个质数，只有当它可以表示成 \(p = 2 ^ \alpha + 1\)，且需要满足多项式的项数 \(< \alpha\) 时才能使用 NTT。\(p\) 后面有个加一，是因为 \(g_n\) 指数的分子上出现了 \(-1\)；\(p - 1\) 需要时二的整数次幂，是因为每次都要除以 \(2\)。

bonus：常用质数及原根

\(p = 998244353 = 119 \times 2 ^ {23} + 1, g = 3\)

\(p = 1004535809 = 479 \times 2 ^ {21} + 1, g = 3\)

void NTT(int *a, int n, int op) {
	if (n == 1) return;
	int a1[n], a2[n];
	rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
	NTT(a1, n >> 1, op), NTT(a2, n >> 1, op);
	int gn = qpow((op == 1 ? g : invg), (mod - 1) / n), gk = 1;
	rop(i, 0, n >> 1) {
		a[i] = (a1[i] + 1ll * gk * a2[i]) % mod;
		a[i + (n >> 1)] = (a1[i] - 1ll * gk * a2[i] % mod + mod) % mod;
		gk = 1ll * gk * gn % mod;
	}
}

NTT 优化：蝴蝶变换

盗大佬一张图

这是进行 NTT 的过程中数组下标的变化。

这样看似乎毫无规律。但是把他们写成二进制：

变换前：\(000 ~ 001 ~ 010 ~ 011 ~ 100 ~ 101 ~ 110 ~ 111\)

变换后：\(000 ~ 100 ~ 010 ~ 110 ~ 001 ~ 101 ~ 011 ~ 111\)

可以发现，就是对每个下标进行了二进制翻转。

因此可以先预处理出每个下标在递归底层对应的新下标。然后从底层往顶层迭代合并。

预处理每个下标的二进制翻转：

rop(i, 0, n) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit;

优化后的 NTT：

void NTT(int *a, int n, int op) {
	rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
		int gn = qpow((op == 1 ? g : invg), (mod - 1) / (mid << 1));
		for (int i = 0, gk = 1; i < n; i += (mid << 1), gk = 1)
			for (int j = 0; j < mid; j ++ , gk = 1ll * gk * gn % mod) {
				int x = a[i + j], y = 1ll * a[i + j + mid] * gk % mod;
				a[i + j] = Mod(x + y), a[i + j + mid] = Mod(x - y);
			}
	}
}

当然了，FFT 也可以用蝴蝶变换来优化。实践证明，速度变快了至少二分之一。

FFT 的迭代实现

void FFT(cp *a, int n, int op) {
	rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
		cp Wn = {cos(pi / mid), op * sin(pi / mid)};
		for (int i = 0; i < n; i += (mid << 1)) {
			cp Wk = {1, 0};
			for (int j = 0; j < mid; j ++ , Wk = Wk * Wn) {
				cp x = a[i + j], y = Wk * a[i + j + mid];
				a[i + j] = x + y, a[i + j + mid] = x - y;
			}
		}
	}
}

任意模数多项式乘法（MTT）

计算 \(f \times g(\bmod ~ p)\) 的结果（\(p\) 不一定能够表示成 \(2 ^ \alpha + 1\) 的形式）。

这个东西有两种做法，但是我只学会了三模 NTT。

首先，卷积之后每个系数最多达到 \(\max \{V\} ^ 2 \times n\) 的大小。拿模板题举例，这个东西就是 \(10 ^ {23}\)。因此只需要找三个模数 \(p_1, p_2, p_3\)，满足 \(p_1p_2p_3 > \max \{V\} ^ 2 \times n\)，然后用他们分别 NTT，最后再 CRT / exCRT 合并即可。

void NTT(int *a, int n, int op, int p) {
	rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
		int gn = qpow((op == 1 ? g : invg), (p - 1) / (mid << 1), p);
		for (int i = 0, gk = 1; i < n; i += (mid << 1), gk = 1)
			for (int j = 0; j < mid; j ++ , gk = gk * gn % p) {
				int x = a[i + j], y = a[i + j + mid] * gk % p;
				a[i + j] = Mod(x + y, p), a[i + j + mid] = Mod(x - y, p);
			}
	}
}
int CRT(int a, int b, int c) {
	int k = Mod(b - a, p2) * inv1 % p2;
	int x = k * p1 + a;
	k = Mod(c - x, p3) * inv2 % p3;
	x = (x + k * p1 % p * p2 % p) % p;
	return x;
}
signed main() {
	read(n, m, p);
	rep(i, 0, n) read(a[i]);
	rep(i, 0, m) read(b[i]);
	memcpy(a1, a, sizeof a); memcpy(a2, a, sizeof a);
	memcpy(a3, a, sizeof a); memcpy(b1, b, sizeof b);
	memcpy(b2, b, sizeof b); memcpy(b3, b, sizeof b);
	m = n + m; n = 1; int bit = 0;
	while (n <= m) n <<= 1, bit ++ ; bit -- ;
	rop(i, 0, n) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit;
	invg = inv(g, p1); NTT(a1, n, 1, p1); NTT(b1, n, 1, p1);
	invg = inv(g, p2); NTT(a2, n, 1, p2); NTT(b2, n, 1, p2);
	invg = inv(g, p3); NTT(a3, n, 1, p3); NTT(b3, n, 1, p3);
	rop(i, 0, n) a1[i] = 1ll * a1[i] * b1[i] % p1;
	rop(i, 0, n) a2[i] = 1ll * a2[i] * b2[i] % p2;
	rop(i, 0, n) a3[i] = 1ll * a3[i] * b3[i] % p3;
	invg = qpow(g, p1 - 2, p1); NTT(a1, n, -1, p1);
	invg = qpow(g, p2 - 2, p2); NTT(a2, n, -1, p2);
	invg = qpow(g, p3 - 2, p3); NTT(a3, n, -1, p3);
	inv1 = inv(n, p1); inv2 = inv(n, p2); inv3 = inv(n, p3);
	rop(i, 0, n) a1[i] = a1[i] * inv1 % p1, a2[i] = a2[i] * inv2 % p2, a3[i] = a3[i] * inv3 % p3;
	inv1 = inv(p1, p2); inv2 = inv(p1 * p2 % p3, p3);
	rep(i, 0, m) printf("%lld ", CRT(a1[i], a2[i], a3[i]));
	return 0;
}

这里选择的模数为：\(p_1 = 998244353, p_2 = 469762049, p_3 = 1004535809\)。他们的原根都为 \(g = 3\)。

多项式求逆

求多项式 \(f(x)\) 的逆元 \(f^{-1}(x)\)。\(f(x)\) 的逆元定义为满足 \(f(x) g(x) = 1(\bmod \ x ^ n)\) 的多项式 \(g(x)\)。

使用倍增法即可求出多项式的逆元。

首先假设已经求出了 \(f(x) g'(x) \equiv 1(\bmod \ x ^ {\lceil \frac{n}{2} \rceil})\)。再假设 \(f(x) \bmod \ x ^ n\) 意义下逆元为 \(g(x)\)。那么有：

\[g(x) \equiv g'(x) (\bmod \ x ^ {\lceil \frac{n}{2} \rceil}) \]

\[(g(x) - g'(x)) \equiv 0 (\bmod \ x ^ {\lceil \frac{n}{2} \rceil}) \]

\[(g(x) - g'(x)) ^ 2 \equiv 0 (\bmod \ x ^ n) \]

\[g^2(x) - 2 g(x) g'(x) + g'^2(x) \equiv 0 (\bmod \ x ^ {n}) \]

两边同时乘以 \(f(x)\)，可以得到：

\[g(x) - 2 g'(x) + f(x) g'^2(x) \equiv 0 (\bmod \ x ^ n) \]

\[g(x) \equiv 2g'(x) - f(x) g'^2(x) (\bmod \ x ^ n) \]

\[g(x) \equiv g'(x) (2 - f(x)g'(x)) (\bmod \ x ^ n) \]

倍增求即可。

void Inv(int *f, int *g, int n) {
	if (n == 1) {
		g[0] = inv(f[0]); return;
	} Inv(f, g, (n + 1) >> 1);
	int m = 1, bit = 0;
	while (m < (n << 1)) m <<= 1, bit ++ ; bit -- ;
	rop(i, 0, n) tmp[i] = f[i];
	rop(i, n, m) tmp[i] = 0;
	rev = new int [m + 5]();
	rop(i, 1, m) rev[i] = (rev[i >> 1] >> 1) | (i & 1) << bit;
	NTT(tmp, m, 1); NTT(g, m, 1);
	rop(i, 0, m) g[i] = (2ll - 1ll * g[i] * tmp[i] % p + p) % p * g[i] % p;
	NTT(g, m, -1); int invn = inv(m);
	rop(i, 0, m) g[i] = 1ll * g[i] * invn % p;
	rop(i, n, m) g[i] = 0; 
	free(rev); rev = NULL;
}