约定：设 $\deg(A(x))$ 为多项式 $A(x)$ 的次数。
FFT 和 NTT
可以参见 Miskcoo’s Blog
多项式求逆
定义： 对于 $n$ 次多项式 $A(x)$ ，求出多项式 $B(x)$ ，满足 $A(x) \times B(x) \equiv 1 \pmod{x^p}$.
用数学归纳法的思路。
在 $ \text{mod } x $ 的意义下， $A(x)$ 的逆多项式为 ${a_0}^{-1}$ . 也就是说，要能多项式求逆，首先常数项必须存在逆。在 $\bmod 998244353$ 意义下这等价于 $a_0 \ne 0$.
假定我们已经知道了 $A(x)$ 在 $\text{mod }x^t$ 意义下的逆 $B_0(x)$ ，想要求它在 $\text{mod }x^{2t}$ 意义下的逆 $B(x)$.
由逆多项式定义， $A(x) \times B(x) \equiv 1 \pmod{x^{2t}}$ ，也就是 $1 \dots x^{2t-1}$ 次项系数都是 $0$ ，只有常数项是 $1$ .
那么显然就有
\[A(x) \times B(x) \equiv 1 \pmod{x^{t}}\]而由定义又有
\[A(x) \times B_0(x) \equiv 1 \pmod{x^{t}}\]两式相减，得
\[A(x) \times [B(x)-B_0(x)] \equiv 0 \pmod{x^{t}}\]显然 $A(x) \ne 0$ ，否则逆多项式不存在。于是一定有
\[B(x) - B_0(x) \equiv 0 \pmod{x^{t}}\]这等价于 $x^t \mid [B(x) - B0(x)]$ .
平方后，有
\[x^{2t} \mid [B(x) - B_0(x)]^2 \Leftrightarrow [B(x) - B_0(x)]^2 \equiv 0 \pmod{x^{2t}} \\ B(x)^2 \equiv 2B(x) \times B_0(x) - B_0(x)^2 \pmod{x^{2t}}\]同乘 $A(x)$ 后： \(B(x) \equiv 2B_0(x) - B_0(x)^2 \times A(x) \pmod{x^{2t}}\)
稍微变形一下，就是： \(A(x) \times B_0(x) \equiv 1 \pmod{ x^{\lceil \frac{p}{2} \rceil}} \\ \Rightarrow B(x) \equiv B_0(x)\times(2 - A(x) \times B_0(x) ) \pmod{x^{p}}\)
倍增即可。（想一想，为什么是上取整）
Code:
void poly_inv(int deg, int A[], int B[]) { // mod x^{deg}，调用前要保证B[]为空
static int tmp[MAXN+10];
if(deg == 1) {
B[0] = mod_inv(A[0]);
} else {
poly_inv((deg + 1) >> 1, A, B);
int len;
for(len = 1; len < 2 * deg; len <<= 1); // 注意：相乘后次数>=n，所以开2倍
copy(A, A+deg, tmp);
fill(tmp+deg, tmp+len, 0);
NTT(len, tmp); NTT(len, B);
for(int i = 0; i < len; i++)
B[i] = mul(B[i], sub(2, mul(tmp[i], B[i])));
NTT(len, B, true);
fill(B+deg, B+len, 0);
}
}
多项式除法 / 多项式取模
定义：对于 $n$ 次多项式 $A(x)$ 和 $m$ 次多项式 $B(x)$ ，求出 $n-m$ 次多项式 $Q(x)$ 和小于 $m$ 次的多项式 $R(x)$，使得它们满足下列关系： \(A(x) = Q(x) \times B(x) + R(x)\)
思路是考虑用模去掉余数的影响。
取 $A(x)$ 的系数翻转 $A^R(x) = x^nA(\frac{1}{x})$ ，于是：
\[\begin{align*} A^R(x) &= x^n \times Q(\frac{1}{x}) \times B(\frac{1}{x}) + x^n \times R(\frac{1}{x}) \\ &= \left(x^{n-m}Q(\frac{1}{x})\right) \times \left(x^mB(\frac{1}{x})\right) + x^n \times R(\frac{1}{x}) \\ &= Q^R(x) \times B^R(x) + x^n \times R(\frac{1}{x}) \end{align*}\]显然 $x^n \times R(1/x)$ 的系数非 $0$ 项次数在 $[n-m+1, n] \cap \mathbb{Z}$ 内。于是放在 $\bmod{x^{n-m+1}}$ 意义下即可消除其影响。
\[\begin{align*} A^R(x) &\equiv Q^R(x) \times B^R(x) &\pmod{x^{n-m+1}} \\ Q^R(x) &\equiv A^R(x) \times \left(B^R(x)\right)^{-1} &\pmod{x^{n-m+1}} \end{align*}\]刚好 $Q^R(x)$ 的次数小于 $n-m+1$ ，故可以完整地求出 $Q^R(x)$ ，进而求出 $Q(x)$，剩下的也好办了。
Code:
void poly_sub(int deg, int A[], int B[], int C[]) {
for(int i = 0; i <= deg; i++) C[i] = sub(A[i], B[i]);
}
void poly_div(int n, int m, int A0[], int B0[], int Q[], int R[]) {
static int A[MAXN+10], B[MAXN+10], tmp[MAXN+10];
CLEAR(A); CLEAR(B); CLEAR(tmp);
copy(A0, A0+n+1, A); copy(B0, B0+m+1, B);
reverse(A, A+n+1); reverse(B, B+m+1);
poly_inv(n-m+1, B, tmp);
conv(n-m, A, tmp, Q);
fill(Q+n-m+1, Q+n+1, 0);
reverse(Q, Q+n-m+1);
copy(A0, A0+n+1, A); copy(B0, B0+m+1, B);
CLEAR(tmp);
conv(n, B, Q, tmp);
poly_sub(n, A, tmp, R);
fill(R+m, R+n+1, 0);// n+1! past-end iterator!
}
多项式牛顿迭代
数域为 $\mathbb{R}$ 的牛顿迭代
用途：求解函数 $f(x)$ 的零点。
我们考虑泰勒展开：
\[g(x) = \sum_{i=0}^{\infty}\frac{f^{(i)}(x_0)}{i!} (x-x_0)^i\]取前 $2$ 项来拟合 $f(x)$，求其零点，即为：
\[g(x) \approx f(x_0) + f'(x_0)\cdot(x-x_0) = 0\Rightarrow x = x_0 - \frac{f(x_0)}{f'(x_0)}\]不断迭代即可接近零点。
多项式环上的牛顿迭代
用途：已知 $\deg(G(x)) = t$. 求解满足下列方程的 $F(x)$：
\[G(F(x)) \equiv 0 \pmod{x^n}\]考虑类似上文的迭代，每次精度翻倍，即增加了一倍的系数（对于高次项）
当 $n = 1$ 时必须有 $G(F(x)) \equiv 0 \pmod{x}$，这是数学归纳法的基础。
假设我们已知 $F_0(x)$ ，且 $\deg F_0(x) = \left\lceil\frac{n}{2}\right\rceil$. 考虑在 $F_0(x)$ 处对 $G(x)$ 进行泰勒展开：
\[0 = G(F_0(x)) + G'(F_0(x)) \cdot (F(x)-F_0(x)) + \frac{G''(F_0(x))}{2} \cdot \left(F(x)-F_0(x)\right) ^2 + \cdots\]注意 $\left(F(x) - F_0(x)\right)^2$ 这一项：由于 $F(x) \equiv F_0(x) \pmod{x^{\lceil n/2 \rceil}}$ ，在 $\bmod{x^n}$ 意义下，$F(x)-F_0(x)$ 系数非 $0$ 项次数属于 $\left[\lceil n/2 \rceil, n\right] \cap \mathbb{Z}$.
于是 $\left(F(x) - F_0(x)\right)^2$ 的非 $0$ 项中，最低的一项次数 $\ge n$. 故它即之后的项对答案没有贡献。
故最终的式子是：
\[F(x) = F_0(x) - \frac{G(F_0(x))}{G'(F_0(x)))}\]多项式 $\ln$
定义：求出多项式 $F(x)$ ，满足 $F(x) \equiv \ln G(x) \pmod{x^n}$.
这里我们不需要采用多项式算法，转而使用微积分解决。
\[\begin{align*} \frac{\text{d}}{\text{d}x} \ln G(x) &= \frac{1}{G(x)} \cdot \frac{\text{d}}{\text{d}x} G(x) \\ \int \text{d}\ln G(x) &= \int \frac{1}{G(x)} \cdot \text{d} G(x) \\ \ln G(x) &= \int \frac{G'(x)}{G(x)} \text{d}x \end{align*}\]其中 $\frac{G’(x)}{G(x)} = G’(x)\times G^{-1}(x)$，用多项式求逆解决。下面的问题就是如何实现多项式求导和积分。
求导：
\[\frac{\text{d}}{\text{d}x}\sum_{i=0}^n a_ix^i = \sum_{i=0}^{n-1}(i+1) \cdot a_{i+1}x^i\]积分：
\[\int \sum_{i=0}^n a_i x^i = \sum_{i=1}^{n+1} \frac{a_{i-1}}{i} x^i\]这些都可以 $O(n)$ 实现。
Code:
void poly_int(int deg, int A[], int B[]) {
B[0] = 0;
for(int i = 1; i <= deg + 1; i++)
B[i] = mul(A[i - 1], mod_inv(i));
}
void poly_d(int deg, int A[], int B[]) {
for(int i = 0; i <= deg - 1; i++)
B[i] = mul(i + 1, A[i + 1]);
B[deg] = 0;
}
void poly_ln(int deg, int A0[], int B[]) { // mod x^{deg}
static int A[MAXN + 10], A1[MAXN + 10], tmp[MAXN + 10];
fill(A, A + MAXN, 0); fill(A1, A1 + MAXN, 0); fill(tmp, tmp + MAXN, 0);
copy(A0, A0+deg, A);
poly_d(deg - 1, A, A1);
poly_inv(deg, A, tmp);
copy(tmp, tmp + deg, A);
conv(deg-1, A, A1, tmp);
poly_int(deg-1, tmp, B);
}
多项式 $\exp$
定义：求出多项式 $F(x)$，满足 $F(x)\equiv e^{G(x)} \pmod{x^n} \pmod{x^n}$.
变形为 $\ln F(x) - G(x) = 0$ ，然后直接套用多项式牛顿迭代，则有：
\[F(x) = F_0(x) + (G(x)-\ln F_0(x)) \times F_0(x)\]使用多项式 $\ln$ 计算 $\ln F_0(x)$ 即可。
Code:
void poly_exp(int deg, int A[], int B[]) { // mod x^{deg}
static int tmp[MAXN + 10];
fill(tmp, tmp + MAXN, 0);
if(deg == 1) {
B[0] = 1;
} else {
int half = (deg + 1) >> 1;
poly_exp(half, A, B);
fill(B+half, B+deg, 0);
poly_ln(deg, B, tmp);
poly_sub(deg - 1, A, tmp, tmp);
conv(deg - 1, tmp, B, tmp);
poly_add(deg - 1, B, tmp, B);
}
}
多项式开方 / 多项式任意次幂
定义：求出多项式 $F(x)$，满足 $F(x) \equiv {G(x)}^k \pmod{x^n}$.
特殊地，定义多项式开方为 $k=\dfrac{1}{2}$，即 $F(x) \equiv \sqrt{G(x)}$ 的情形。
显然有 $F(x) \equiv \exp(k\ln(G(x)))$，于是应用多项式 $\ln$ 和多项式 $\exp$ 即可解决。
多项式模板大全
#include <bits/stdc++.h>
using namespace std;
namespace polynomial {
const int MAXN = 100000 << 2, MOD = 998244353, G = 3;
int p_omega[MAXN + 10], p_inv[MAXN + 10];
inline int pls(int a, int b) {
int ret = (a + b);
if(ret >= MOD) ret -= MOD;
return ret;
}
inline int dec(int a, int b) {
int ret = (a - b);
if(ret < 0) ret += MOD;
return ret;
}
inline int mul(int a, int b) { return (1ll * a) * b % MOD; }
inline int fpow(int a, int x) {
int ret = 1;
while(x) {
if(x & 1) ret = mul(ret, a);
x >>= 1;
a = mul(a, a);
}
return ret;
}
inline int mod_inv(int x) { return fpow(x, MOD - 2); }
//-------------------NTT----------------------
void init_omega() {
for(int step = 1; step <= MAXN; step <<= 1) {
p_omega[step] = fpow(G, (MOD - 1) / step);
p_inv[step] = fpow(G, (MOD - 1) / step * (step - 1));
}
}
struct _omega_init_t {
_omega_init_t() { init_omega(); }
} _omega_init;
inline int get_rt(int step, bool idft) {
return !idft ? p_omega[step] : p_inv[step];
}
inline void NTT(int len, int A[], bool idft = false) {
static int rev[MAXN + 10];
for(int i = 0; i < len; i++)
rev[i] = (((rev[i >> 1] >> 1) | (len >> (i & 1))) & (len - 1));
for(int i = 0; i < len; i++)
if(rev[i] > i) swap(A[i], A[rev[i]]);
for(int step = 1; step < len; step <<= 1) {
int rt = get_rt(step << 1, idft);
for(int i = 0; i < len; i += (step << 1)) {
for(int j = 0, omega = 1; j < step;
j++, omega = mul(omega, rt)) {
int t = mul(A[i + j + step], omega);
A[i + j + step] = dec(A[i + j], t);
A[i + j] = pls(A[i + j], t);
}
}
}
if(idft) {
int t = mod_inv(len);
for(int i = 0; i < len; i++) A[i] = mul(A[i], t);
}
}
inline void conv(int deg, int A0[], int B0[], int C[]) {
static int A[MAXN + 10], B[MAXN + 10];
int len;
copy(A0, A0+deg+1, A), copy(B0, B0+deg+1, B);
for(len = 1; len < deg * 2 + 1; len <<= 1);
fill(A+deg+1, A+len, 0), fill(B+deg+1, B+len, 0); // !!
NTT(len, A), NTT(len, B);
for(int i = 0; i < len; i++) C[i] = mul(A[i], B[i]);
NTT(len, C, true);
}
//-------------------Polynomial Add & Substract----------
void poly_add(int deg, int A[], int B[], int C[]) {
for(int i = 0; i <= deg; i++) C[i] = pls(A[i], B[i]);
}
void poly_sub(int deg, int A[], int B[], int C[]) {
for(int i = 0; i <= deg; i++) C[i] = dec(A[i], B[i]);
}
//------------Polynomial Integral & Derivative-----------
void poly_int(int deg, int A[], int B[]) {
B[0] = 0;
for(int i = 1; i <= deg + 1; i++)
B[i] = mul(A[i - 1], mod_inv(i));
}
void poly_d(int deg, int A[], int B[]) {
for(int i = 0; i <= deg - 1; i++)
B[i] = mul(i + 1, A[i + 1]);
B[deg] = 0;
}
//-------------------Polynomial inverse------------------
void poly_inv(int deg, int A[], int B[]) { // mod x^{deg}
int tmp[MAXN + 10];
if(deg == 1) {
assert(A[0] != 0);
B[0] = mod_inv(A[0]);
} else {
poly_inv((deg + 1) >> 1, A, B);
int len;
for(len = 1; len < 2 * deg; len <<= 1); // !!
copy(A, A+deg, tmp); fill(tmp+deg, tmp+len, 0);
NTT(len, tmp); NTT(len, B);
for(int i = 0; i < len; i++)
B[i] = mul(B[i], dec(2, mul(tmp[i], B[i])));
NTT(len, B, true);
fill(B+deg, B+len, 0); // !!
}
}
//--------------------Polynomial division-----------------
void poly_div(int n, int m, int A0[], int B0[], int Q[], int R[]) { // mod x^{deg}
int A[MAXN + 10], B[MAXN + 10], tmp[MAXN + 10];
fill(A, A+MAXN, 0); fill(B, B+MAXN, 0); fill(tmp, tmp+MAXN, 0);
copy(A0, A0+n+1, A), copy(B0, B0+m+1, B);
reverse(A, A+n+1), reverse(B, B+m+1);
poly_inv(n-m+1, B, tmp);
conv(n-m, A, tmp, Q);
fill(Q+n-m+1, Q+n+1, 0);
reverse(Q, Q+n-m+1);
copy(A0, A0+n+1, A), copy(B0, B0+m+1, B);
conv(n, B, Q, tmp);
poly_sub(n, A, tmp, R);
fill(R+m, R+n+1, 0);
}
//--------------------Polynomial ln------------------------
void poly_ln(int deg, int A0[], int B[]) { // mod x^{deg}
static int A[MAXN + 10], A1[MAXN + 10], tmp[MAXN + 10];
fill(A, A + MAXN, 0); fill(A1, A1 + MAXN, 0); fill(tmp, tmp + MAXN, 0);
copy(A0, A0+deg, A);
poly_d(deg - 1, A, A1);
poly_inv(deg, A, tmp);
copy(tmp, tmp + deg, A);
conv(deg-1, A, A1, tmp);
poly_int(deg-1, tmp, B);
}
//--------------------Polynomial exp-----------------------
void poly_exp(int deg, int A[], int B[]) { // mod x^{deg}
static int tmp[MAXN + 10];
fill(tmp, tmp + MAXN, 0);
if(deg == 1) {
B[0] = 1;
} else {
int half = (deg + 1) >> 1;
poly_exp(half, A, B);
fill(B+half, B+deg, 0);
poly_ln(deg, B, tmp);
poly_sub(deg - 1, A, tmp, tmp);
conv(deg - 1, tmp, B, tmp);
poly_add(deg - 1, B, tmp, B);
}
}
}
