多项式全家桶
Galx_Trail
·
·
算法·理论
前言：lyy 还是太强了 /bx
前置芝士
代数基本定理
- 一个 d 次多项式可以被 d + 1 个点唯一确定。
- 还有要用的后面补。
虚数，复数杂烩
基础定义：i^2 = -1。
复平面：在一维实数轴上扩展一个虚数轴。
复数的表示
- 代数法：a + bi，其中 a 为实部，bi 为虚部。
- 极角表示：复平面上一个点到原点 O 的向量，其模长为 r，辐角为 \theta，则可以表示为：r \times \cos\theta + r \times i\sin\theta。
在复数计算上，加法相当于向量相加，乘法相当于模长相乘、辐角相加，这在后面理解 FFT 时非常关键。
欧拉公式
e^{i\theta} = \cos\theta + i\sin\theta
从复平面的视角理解，其几何意义就是当 \theta \in [0, 2\pi] 时，e^{i\theta} 刚好画出一个单位圆。
单位复数根
求解方程：
z^n = 1
在实数域上只有 1（偶次方还有 -1）。根据代数基本定理，在复数域上有 n 个根。
$$z^n = e^{i2k\pi} \implies z = e^{i\frac{2k\pi}{n}}$$
这 $n$ 个根均匀分布在单位圆上，令逆时针方向从 $(1,0)$ 开始的第一个点为主 $n$ 次单位根，记为 $\omega_n$。剩下的所有根都是 $\omega_n$ 多次旋转相同角度得到的，即：$\omega_n^0, \omega_n^1, \omega_n^2, \dots, \omega_n^{n-1}$。
#### 一些引理
**1. 消去引理**
$$\omega_{dn}^{dk} = \omega_n^k$$
复平面上对应的旋转角度 $\theta$ 不变，因此数值相等。
**2. 折半引理**
如果 $n$ 是偶数，那么这 $n$ 个单位根的平方，恰好对应 $\frac{n}{2}$ 个 $\frac{n}{2}$ 次单位根，且每个出现 $2$ 次。
平方意味着旋转角度翻倍。以 $n = 8$ 为例，上下部分各 $4$ 个点，平方后上面 $4$ 个点与下面 $4$ 个点旋转后的位置完全重合。这是 FFT 能做到 $O(n\log n)$ 的核心分治依据。
**3. 求和引理**
$$\sum_{j=0}^{n-1} (\omega_n^k)^j = 0 \quad (k \not\equiv 0 \pmod n)$$
从复平面视角看，根据 $n$ 次单位根的定义，$n$ 个根构成了一个中心在原点的正 $n$ 边形顶点。由于其具备旋转对称性，所有向量和一定为 $0$。
---
## 拉格朗日插值
### 一般公式
已知 $m$ 个点 $(x_1, y_1), (x_2, y_2), \dots, (x_m, y_m)$（其中 $x_i$ 互不相同），寻找一个最高次数为 $m - 1$ 次的多项式 $f(x)$ 使得：$\forall i, f(x_i) = y_i$。
构造形式如下：
$$\sum_{i=1}^{m} y_i \times \prod_{j \not= i} \frac{x - x_j}{x_i - x_j}$$
当代入 $x = x_i$ 时，除第 $i$ 项外的所有项分子均包含 $(x_i - x_i) = 0$，而第 $i$ 项的分式为 $1$，整体结果恰好为 $y_i$。
#### 模板实现（[P4781](https://www.luogu.com.cn/problem/P4781)）
```cpp
constexpr ll mod = 998244353;
ll fpow(ll a, ll b) {
ll res = 1;
a %= mod;
while (b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
ll inv(ll x) {
return fpow(x, mod - 2);
}
void Main() {
ll k;
int n; cin >> n >> k;
vec<pair<ll, ll>> a(n + 1);
for (int i = 1; i <= n; ++i) cin >> a[i].first >> a[i].second;
ll ans = 0;
for (int i = 1; i <= n; ++i) {
ll s1 = a[i].second % mod, s2 = 1;
for (int j = 1; j <= n; ++j) {
if (i == j) continue;
s1 = (s1 * ((k - a[j].first) % mod + mod)) % mod;
s2 = (s2 * ((a[i].first - a[j].first) % mod + mod)) % mod;
}
ans = (ans + s1 * inv(s2) % mod) % mod;
}
cout << ans;
}
```
### 拓展问题
[CF622F](https://codeforces.com/problemset/problem/622/F)：定义函数 $S_k(n) = \sum_{i=1}^{n} i^k$，求 $S_k(n) \pmod{10^9+7}$（其中 $n \le 10^{18}, k \le 10^6$）。
设 $f(n) = S_k(n)$，前向差分 $\Delta f(n) = f(n + 1) - f(n) = (n + 1)^k$。由于 $\Delta f(n)$ 为 $k$ 次多项式，根据有限差分性质，$S_k(n)$ 是一个 $k + 1$ 次多项式。
确定该多项式需要 $k + 2$ 个点。当点值横坐标连续时（取 $x_i = i$），可将复杂度优化至 $O(k)$：
$$\sum_{i=1}^{k+2} y_i \times \prod_{j \not= i} \frac{n - j}{i - j}$$
* **分子部分**：$\frac{\prod_{j=1}^{k+2} (n - j)}{n - i}$，可通过前缀积和后缀积在 $O(k)$ 内预处理。
* **分母部分**：
$$\prod_{j=1, j\not=i}^{k+2} (i - j) = (i - 1)! \times (-1)^{k+2-i} \times (k + 2 - i)!$$
使用阶乘及阶乘逆元即可 $O(1)$ 求解单项分母。
```cpp
constexpr ll mod = 1e9 + 7, N = 1e6 + 5;
ll fpow(ll a, ll b) {
ll res = 1;
a %= mod;
while (b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int primes[N], cnt;
ll pk[N];
bool is[N];
void init() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout.tie(nullptr);
pk[1] = 1;
}
void Main() {
int n, k;
cin >> n >> k;
for (int i = 2; i < N; ++i) {
if (!is[i]) {
primes[++cnt] = i;
pk[i] = fpow(i, k);
}
for (int j = 1; j <= cnt && i * primes[j] < N; ++j) {
is[i * primes[j]] = true;
pk[i * primes[j]] = pk[i] * pk[primes[j]] % mod;
if (i % primes[j] == 0) break;
}
}
vec<ll> y(k + 5, 0);
for (int i = 1; i < k + 5; ++i) {
y[i] = (y[i - 1] + pk[i]) % mod;
}
if (n <= k + 2) {
cout << y[n];
return;
}
int m = k + 2;
vec<ll> pre(m + 3, 0), suf(m + 3, 0);
pre[0] = 1, suf[m + 1] = 1;
for (int i = 1; i <= m; ++i) {
pre[i] = pre[i - 1] * (n - i) % mod;
}
for (int i = m; i >= 1; --i) {
suf[i] = suf[i + 1] * (n - i) % mod;
}
vec<ll> inv(m + 3, 0), fact(m + 3, 0);
fact[0] = 1;
for (int i = 1; i <= m; ++i) fact[i] = fact[i - 1] * i % mod;
inv[m] = fpow(fact[m], mod - 2);
for (int i = m - 1; i >= 1; --i) inv[i] = inv[i + 1] * (i + 1) % mod;
inv[0] = 1;
ll ans = 0;
for (int i = 1; i <= m; ++i) {
ll sum1 = pre[i - 1] * suf[i + 1] % mod;
ll sum2 = inv[i - 1] * inv[m - i] % mod;
ll sum3 = y[i] * sum1 % mod * sum2 % mod;
if ((m - i) % 2 == 0) {
ans = (ans + sum3) % mod;
} else {
ans = ((ans - sum3) % mod + mod) % mod;
}
}
cout << ans;
}
```
#### [P5364 [SNOI2017] 礼物](https://www.luogu.com.cn/problem/P5364)
由题意：
$$f_i = f_{i-1} + i^k, \quad s_i = s_{i-1} + f_i = 2s_{i-1} + i^k$$
展开后项：
$$i^k = ((i - 1) + 1)^k = \sum_{j=0}^{k} \binom{k}{j} (i - 1)^j$$
构造状态向量：$[s_i, i^0, i^1, \dots, i^k]^T$。转移矩阵大小为 $(k + 2) \times (k + 2)$，使用矩阵快速幂可在 $O(k^3 \log n)$ 完成递推。
```cpp
constexpr ll mod = 1e9 + 7;
int sz;
ll fpow(ll a, ll b) {
ll res = 1;
a %= mod;
while (b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
struct Mat {
ll mat[15][15];
Mat() {
for (int i = 0; i <= sz; ++i)
for (int j = 0; j <= sz; ++j) mat[i][j] = 0;
}
};
Mat mul(const Mat &a, const Mat &b) {
Mat c;
for (int i = 0; i < sz; ++i) {
for (int kk = 0; kk < sz; ++kk) {
if (!a.mat[i][kk]) continue;
for (int j = 0; j < sz; ++j) {
c.mat[i][j] = (c.mat[i][j] + a.mat[i][kk] * b.mat[kk][j] % mod) % mod;
}
}
}
return c;
}
Mat Mat_fpow(Mat a, ll b) {
Mat res;
for (int i = 0; i < sz; ++i) res.mat[i][i] = 1;
while (b) {
if (b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}
void Main() {
ll n; int k;
cin >> n >> k;
sz = k + 2;
vec<vec<ll>> C(15, vec<ll>(15, 0));
for (int i = 0; i < 15; ++i) {
C[i][0] = 1;
for (int j = 1; j <= i; ++j) {
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
}
}
Mat base;
base.mat[0][0] = 2;
for (int j = 0; j <= k; ++j) base.mat[0][j + 1] = C[k][j];
for (int i = 0; i <= k; ++i) {
for (int j = 0; j <= i; ++j) {
base.mat[i + 1][j + 1] = C[i][j];
}
}
Mat res = Mat_fpow(base, n - 1);
ll sum = res.mat[0][1];
cout << (sum + fpow(n, k)) % mod;
}
```
---
## FFT
### 核心思想
计算卷积 $C = A * B$ 时，暴力系数相乘复杂度为 $O(n^2)$。FFT 将多项式转化为点值表达式，进行 $O(n)$ 点值乘法后再插值还原。
设多项式 $A(x) = \sum_{j=0}^{n-1} a_j x^j$，将其按奇偶项拆分：
$$A^{[0]}(x) = a_0 + a_2x + a_4x^2 + \dots$$
$$A^{[1]}(x) = a_1 + a_3x + a_5x^2 + \dots$$
$$A(x) = A^{[0]}(x^2) + x A^{[1]}(x^2)$$
代入 $\omega_n^k$ 与 $\omega_n^{k+n/2}$，由折半引理得：
$$y_k = A(\omega_n^k) = A^{[0]}(\omega_{n/2}^k) + \omega_n^k A^{[1]}(\omega_{n/2}^k)$$
$$y_{k+n/2} = A(\omega_n^{k+n/2}) = A^{[0]}(\omega_{n/2}^k) - \omega_n^k A^{[1]}(\omega_{n/2}^k)$$
这就是蝴蝶操作。
### 逆变换（IFFT）
变换矩阵为范德蒙德矩阵 $V_n$，其逆矩阵元素满足 $(V_n)^{-1}_{j,k} = \frac{1}{n} \omega_n^{-jk}$。因此 IFFT 只需将单位根替换为其共轭复数 $\omega_n^{-k}$，并在最终结果上除以 $n$ 即可。
### 迭代实现（雷德算法 / 位逆序置换）
通过二进制位反转预处理下标，自底向上合并子问题，可避免递归常数与爆栈风险。
```cpp
const lb pi = acos(-1.0);
constexpr int N = 4e6 + 5;
struct Complex {
lb r, i;
Complex(lb rr = 0, lb ii = 0) : r(rr), i(ii) {}
Complex operator+(const Complex &other) const { return {r + other.r, i + other.i}; }
Complex operator-(const Complex &other) const { return {r - other.r, i - other.i}; }
Complex operator*(const Complex &other) const {
return {r * other.r - i * other.i, r * other.i + i * other.r};
}
};
Complex a[N], b[N];
int rev[N];
void FFT(Complex *A, int lim, int op) {
for (int i = 0; i < lim; ++i) {
if (i < rev[i]) swap(A[i], A[rev[i]]);
}
for (int m = 1; m < lim; m <<= 1) {
Complex wn(cos(pi / m), op * sin(pi / m));
for (int j = 0; j < lim; j += (m << 1)) {
Complex w(1, 0);
for (int k = 0; k < m; ++k) {
Complex x = A[j + k], y = A[j + k + m] * w;
A[j + k] = x + y;
A[j + k + m] = x - y;
w = w * wn;
}
}
}
if (op == -1) {
for (int i = 0; i < lim; ++i) A[i].r /= lim;
}
}
void Main() {
int n, m;
cin >> n >> m;
for (int i = 0; i <= n; ++i) cin >> a[i].r;
for (int j = 0; j <= m; ++j) cin >> b[j].r;
int lim = 1, len = 0;
int sz = n + m + 1;
while (lim < sz) {
lim <<= 1;
++len;
}
for (int i = 0; i < lim; ++i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
}
FFT(a, lim, 1);
FFT(b, lim, 1);
for (int i = 0; i < lim; ++i) a[i] = a[i] * b[i];
FFT(a, lim, -1);
for (int i = 0; i < sz; ++i) cout << (int)round(a[i].r) << ' ';
}
```
---
## NTT
NTT 使用模 $p$ 的原根 $g$ 替代单位复数根 $\omega_n$，构造 $\omega_n \equiv g^{\frac{p-1}{n}} \pmod p$，在模意义下完全消除浮点数精度误差。
### [分治 NTT（P4721）](https://www.luogu.com.cn/problem/P4721)
求解卷积递推式 $f_i = \sum_{j=1}^i f_{i-j} g_j$。使用 CDQ 分治，先递归计算左区间 $[l, mid]$ 的 $f$，再计算左区间对右区间 $[mid + 1, r]$ 的 NTT 卷积贡献。
```cpp
constexpr ll mod = 998244353, G = 3, GI = (mod + 1) / 3, N = 4e5 + 5;
ll fpow(ll a, ll b) {
ll res = 1;
a %= mod;
while (b) {
if (b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int rev[N];
ll f[N], g[N], a[N], b[N];
void NTT(ll *A, int n, int op) {
for (int i = 0; i < n; ++i)
if (i < rev[i]) swap(A[i], A[rev[i]]);
for (int m = 1; m < n; m <<= 1) {
ll wn = fpow(op == 1 ? G : GI, (mod - 1) / (m << 1));
for (int j = 0; j < n; j += m << 1) {
ll w = 1;
for (int k = 0; k < m; ++k) {
ll x = A[j + k], y = w * A[j + k + m] % mod;
A[j + k] = (x + y) % mod;
A[j + k + m] = ((x - y) % mod + mod) % mod;
w = w * wn % mod;
}
}
}
if (op == -1) {
ll inv = fpow(n, mod - 2);
for (int i = 0; i < n; ++i) A[i] = A[i] * inv % mod;
}
}
void cdq(int l, int r) {
if (l == r) return;
int mid = (l + r) >> 1;
cdq(l, mid);
int lim = 1, len = 0;
while (lim <= (r - l + 1) + (mid - l + 1)) {
lim <<= 1;
++len;
}
for (int i = 0; i < lim; ++i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
a[i] = (i <= mid - l) ? f[l + i] : 0;
b[i] = (i <= r - l) ? g[i] : 0;
}
NTT(a, lim, 1);
NTT(b, lim, 1);
for (int i = 0; i < lim; ++i) a[i] = a[i] * b[i] % mod;
NTT(a, lim, -1);
for (int i = mid + 1; i <= r; ++i) {
f[i] = (f[i] + a[i - l]) % mod;
}
cdq(mid + 1, r);
}
```
### [ABC352G - Socks 3](https://atcoder.jp/contests/abc352/tasks/abc352_g)
计算抽取袜子期望：$E[X] = \sum_{k=0}^{\infty} P(X > k)$。
设总袜子数 $S = \sum A_i$。取 $k$ 只袜子各不相同的概率为：
$$P(X > k) = \frac{[x^k] \prod_{i=1}^n (1 + A_i x)}{\binom{S}{k}}$$
使用分治 NTT 在 $O(n \log^2 n)$ 下计算多项式 $\prod_{i=1}^n (1 + A_i x)$ 展开后的各项系数。
```cpp
namespace Poly {
vec<int> rev;
void NTT(vec<ll> &A, int lim, int op) {
if ((int)rev.size() != lim) {
rev.resize(lim);
int len = __builtin_ctz(lim);
for (int i = 0; i < lim; ++i)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
}
for (int i = 0; i < lim; ++i)
if (i < rev[i]) swap(A[i], A[rev[i]]);
for (int m = 1; m < lim; m <<= 1) {
ll wn = fpow(op == -1 ? GI : G, (mod - 1) / (m << 1));
for (int j = 0; j < lim; j += (m << 1)) {
ll w = 1;
for (int k = 0; k < m; ++k) {
ll x = A[j + k] % mod, y = A[j + k + m] * w % mod;
A[j + k] = (x + y) % mod;
A[j + k + m] = ((x - y) % mod + mod) % mod;
w = w * wn % mod;
}
}
}
if (op == -1) {
ll inv = fpow(lim, mod - 2);
for (int i = 0; i < lim; ++i) A[i] = A[i] * inv % mod;
}
}
vec<ll> multiply(vec<ll> A, vec<ll> B) {
int n = A.size(), m = B.size(), sz = n + m - 1;
int lim = 1;
while (lim <= sz) lim <<= 1;
A.resize(lim, 0); B.resize(lim, 0);
NTT(A, lim, 1); NTT(B, lim, 1);
for (int i = 0; i < lim; ++i) A[i] = A[i] * B[i] % mod;
NTT(A, lim, -1);
A.resize(sz);
return A;
}
} // namespace Poly
```
---
## FWT
FWT 用于解决位运算卷积 $C_k = \sum_{i \oplus j = k} A_i B_j$。由于每一位二进制独立，可按高维前缀和思想分治处理。
以 XOR 为例：
* 变换：$A'_0 = A_0 + A_1, \quad A'_1 = A_0 - A_1
- 逆变换：A_0 = \frac{A'_0 + A'_1}{2}, \quad A_1 = \frac{A'_0 - A'_1}{2}
constexpr ll mod = 998244353, inv2 = (mod + 1) / 2;
void FWT(vec<ll> &A, int lim, int op, int type) {
for (int m = 1; m < lim; m <<= 1) {
for (int j = 0; j < lim; j += (m << 1)) {
for (int k = 0; k < m; ++k) {
ll x = A[j + k], y = A[j + k + m];
if (type == 1) { // OR
A[j + k + m] = (op == 1) ? (y + x) % mod : (y - x + mod) % mod;
} else if (type == 2) { // AND
A[j + k] = (op == 1) ? (x + y) % mod : (x - y + mod) % mod;
} else { // XOR
if (op == 1) {
A[j + k] = (x + y) % mod;
A[j + k + m] = (x - y + mod) % mod;
} else {
A[j + k] = (x + y) % mod * inv2 % mod;
A[j + k + m] = (x - y + mod) % mod * inv2 % mod;
}
}
}
}
}
}
经典应用（CF662C）
给定 n \times m 的 01 矩阵（n \le 20, m \le 10^5），可翻转任意行或列，求矩阵中 1 的最少数量。
- 设行翻转掩码为 mask，某一列的原状态为 x，翻转后变为 x \oplus mask。
- 贡献为 \min(\text{popcount}(x \oplus mask), n - \text{popcount}(x \oplus mask))。
- 令 A_x 为原矩阵中状态 x 的列数，B_y = \min(\text{popcount}(y), n - \text{popcount}(y))。
- 则答案向量为 C_{mask} = \sum_{x \oplus y = mask} A_x B_y，直接套用异或 FWT 即可在 O(n 2^n) 完成计算。