多项式初等函数
前置知识： 多项式基本操作（求导、积分、FFT）和一些数学知识。
Reference
多项式牛顿迭代
给定多项式 \(g(x)\)，已知有多项式 \(f(x)\) 满足 \(g(f(x)) \equiv 0 \pmod{x^n}\)，求出模 \(x^n\) 意义下的 \(x^n\)。
假设我们已经求得了模 \(x^{\lceil \frac{n}{2} \rceil}\) 意义下的解 \(f_0(x)\)，那么将 \(g(f(x))\) 在 \(f_0(x)\) 处Taylor展开得
由于 \(f(x)-f_0(x)\) 的最低非零项的次数至少是 \(\lceil \frac{n}{2} \rceil\)，因此 \((f(x)-f_0(x))^i \bmod x^n\) 只有在 \(i=0,1\) 的时候不为 \(0\)，上式改写为
整理即可得到
预先求出模 \(x^1\) 意义下的 \(f(x)\) 后倍增即可。
多项式初等函数
多项式求逆
给定多项式 \(f(x)\)，求多项式 \(g(x)\) 使得 \(f(x)g(x)\equiv 1 \pmod{x^n}\)，记 \(g(x)=f^{-1}(x)\)。
设方程 \(h(g(x)) \equiv \frac{1}{g(x)}-f(x) \pmod{x^n}\)，先求出 \(\Big([x^0]f(x)\Big)^{-1}\) 作为 \(n=1\) 时的解，然后用牛顿迭代倍增。
具体的，根据多项式牛顿迭代，得到
其中 \(h(g_0(x)) \equiv 0 \pmod{x^{\lceil \frac{n}{2} \rceil}}\)
整理得
倍增求即可，时间复杂度 \(T(n)=T(\frac{n}{2})+\mathcal{O}(n \log n)=\mathcal{O}(n \log n)\)。
注意常数因子对效率的影响：处理倍增式子的时候可以不暴力乘，而是在点值形式下做完运算后在逆变换回去。
多项式开方
给定多项式 \(f(x)\)，求出多项式 \(g(x)\) 使得 \(g^2(x) \equiv f(x) \pmod{x^n}\)。
首先 \(\Big([x^0]g(x)\Big)^2\equiv [x^0]f(x)\)，若 \([x^0]f(x)\) 没有平方根，则无解，否则任意选取一个即可（选取不同的平方根会得到不同的结果）。
假设求出了模 \(x^{\lceil \frac{n}{2} \rceil}\) 意义下的解 \(g_0(x)\)，应用牛顿迭代，设 \(h(g(x)) \equiv g^2(x) - f(x)\)，得
整理得到
时间复杂度 \(T(n)=T\left(\frac{n}{2}\right)+\mathcal{O}(n \log n)=\mathcal{O}(n \log n)\)。
不保证 \(A[0]\) 是一个完全平方数时需要算二次剩余。
固定模数可以不用 Cipolla 解二次函数，只需要用BSGS 求出来一个 \(t\) 使得 \(g^t \equiv k\)，其中 \(g\) 是原根，模数为 \(998244353\) 时 \(g=3\)。二次剩余就是 \(g^{\frac{t}{2}}\)，复杂度是 \(\mathcal{O}(\sqrt{P})\) ，其中 \(P\) 是模数。
多项式除法
给定 \(n\) 次多项式 \(f(x)\) 和 \(m\) 次 \(g(x)\)，求多项式 \(Q(x)\)，\(R(x)\) 使得：
- \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 次数小于 \(m\)。
- \(f(x)=g(x)Q(x)+R(x)\)。
对于 \(n\) 次多项式 \(f(x)\)，设 \(\mathscr{R}f(x)\) 表示反转 \(f(x)\) 的各项系数得到的多项式，即 \(\mathscr{R}f(x)=x^nf(x^{-1})\)。
变一下式子：
考虑在模 \(x^{n-m+1}\) 意义下求出 \(Q(x)\)，然后 \(R(x)=f(x)-g(x)Q(x)\) 就可以得到 \(R(x)\) 了。
发现在模 \(x^{n-m+1}\) 意义下，\(\mathscr{R}R(x)\) 那一项变成了 \(0\)，于是
求出 \(\mathscr{R}g(x)\) 的逆元即可得到 \(\mathscr{R}Q(x)\)，做一次反转变换即可得到 \(Q(x)\)。
时间复杂度 \(\mathcal{O}(n \log n)\)。
多项式 \(\ln\)
对于 \(n\) 次多项式 \(f(x)\)，在模 \(x^n\) 意义下求出 \(\ln f(x)\)，保证 \([x^0]f(x)=1\)。
对 \(\ln f(x)\) 求导，得到
再积分，得到
求导和积分都是线性的，只需要做一个多项式求逆，时间复杂度 \(\mathcal{O}(n \log n)\)。
当 \([x^0]f(x) \ne 1\) 时，\(\ln f(x)\) 不存在，原因是此时取 \(\ln\) 后不再收敛，在模意义下无意义。
多项式 \(\exp\)
给定 \(n\) 次多项式 \(f(x)\)，若 \(\exp f(x)\) 存在，必须满足 \([x^0]f(x)=0\)，否则常数项不收敛。
牛顿，设方程 \(h(g(x)) \equiv \ln g(x)-f(x) \pmod{x^n}\)，得到
整理得
时间复杂度 \(T(n)=T(\frac{n}{2})+\mathcal{O}(n \log n)=\mathcal{O}(n \log n)\)。
多项式指数幂
给定次多项式 \(f(x)\) 和自然数 \(k\)，计算 \(f^k(x) \pmod{x^n}\)。
当 \([x^0]f(x)=1\) 时，可以先 \(\ln\) 再 \(\exp\)，答案就是 \(\exp(k\ln f(x))\)。
当 \([x^0]f(x) \ne 1\) 时，设最低次项为 \(f_ix^i\)，那么
时间复杂度 \(\mathcal{O}(n \log n)\)。
多项式三角函数
根据Euler公式，得到
以及
在模 \(998244353\) 下，\({\rm i}\) 可以取 \(998244352\) 的二次剩余 \(86583718\) 或 \(911660635\)，然后用多项式 \(\exp\) 即可，时间复杂度 \(\mathcal{O}(n \log n)\)。
多项式求逆即可算 \(\tan\)。
多项式反三角函数
求导再积分，有
时间复杂度 \(\mathcal{O}(n \log n)\)。
多项式全家桶
#include <bits/stdc++.h>
using namespace std;
typedef vector<int> poly;
const int N = 5e6 + 10, MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
long long r = 1;
for(; b; b >>= 1, a = a * a % MOD)
if(b & 1) r = r * a % MOD;
return r;
}
namespace My_poly {
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
int rev[N], iv[N];
inline void make_iv() {
// 处理 [1, N - 1] 的逆元并存到 iv[] 中
iv[1] = 1;
for(int i = 2; i < N; i ++)
iv[i] = Minus(0, 1ll * iv[MOD % i] * (MOD / i) % MOD);
}
inline void make_rev(int n) {
for(int i = 1; i < n; i ++)
rev[i] = ((rev[i >> 1] >> 1) | (i & 1) * (n >> 1));
}
inline void NTT(poly &A, int type) {
// 调用前请确保 rev[] 数组的正确
// type == 1:= DFT
// type == -1:= IDFT
static const int g = 3; // 原根
int n = A.size();
for(int i = 0; i < n; i ++)
if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int h = 2; h <= n; h <<= 1) {
long long step = ksm(g, (MOD - 1) / h);
if(type == -1) step = ksm(step, MOD - 2);
for(int i = 0; i < n; i += h)
for(int j = i, mul = 1; j < i + (h >> 1); j ++, mul = mul * step % MOD) {
int u = A[j], v = 1ll * A[j + (h >> 1)] * mul % MOD;
A[j] = Plus(u, v), A[j + (h >> 1)] = Minus(u, v);
}
}
if(type == -1) {
long long mul = ksm(n, MOD - 2);
for(int i = 0; i < n; i ++)
A[i] = A[i] * mul % MOD;
}
}
poly operator + (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Plus(res[i], rhs[i]);
return res;
}
poly operator - (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Minus(res[i], rhs[i]);
return res;
}
poly operator * (const int lhs, const poly &rhs) {
poly res = rhs;
for(int i = 0; i < rhs.size(); i ++)
res[i] = 1ll * res[i] * lhs % MOD;
return res;
}
poly operator * (const poly &lhs, const int rhs) {
poly res = lhs;
for(int i = 0; i < lhs.size(); i ++)
res[i] = 1ll * res[i] * rhs % MOD;
return res;
}
poly operator * (poly A, poly B) {
int h = 1;
while(h <= A.size() + B.size()) h <<= 1;
A.resize(h, 0), B.resize(h, 0);
make_rev(h);
NTT(A, 1), NTT(B, 1);
for(int i = 0; i < h; i ++) A[i] = 1ll * A[i] * B[i] % MOD;
NTT(A, -1); return A;
}
inline poly derivative(poly A) {
// 多项式求导
for(int i = 1; i < A.size(); i ++)
A[i - 1] = 1ll * i * A[i] % MOD;
if(!A.empty()) A.pop_back();
return A;
}
inline poly integrate(poly A) {
// 多项式积分
// 使用前请保证调用过make_iv函数或init函数
A.emplace_back(0);
for(int i = A.size() - 1; i >= 1; i --)
A[i] = 1ll * A[i - 1] * iv[i] % MOD;
A[0] = 0; // 不定积分的常数 C
return A;
}
inline poly inv(poly A, int n) {
// 多项式求逆
// 返回模 x^n 意义下 A 的逆元
int h = 1; while(h < n) h <<= 1; A.resize(h, 0);
if(A.empty()) return poly();
poly res(1, ksm(A[0], MOD - 2));
for(int i = 2; i <= h; i <<= 1) {
poly q(A.begin(), A.begin() + i);
res.resize(2 * i, 0), q.resize(2 * i, 0);
make_rev(2 * i), NTT(res, 1), NTT(q, 1);
for(int j = 0; j < 2 * i; j ++) res[j] = 1ll * res[j] * Minus(2, 1ll * res[j] * q[j] % MOD) % MOD;
NTT(res, -1); res.resize(i);
}
res.resize(n); return res;
}
inline poly ln(poly A, int n) {
// 多项式 ln，返回 ln A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 1
A.resize(n, 0);
A = derivative(A) * inv(A, n); A.resize(n);
return integrate(A);
}
inline poly exp(poly A, int n) {
// 多项式 exp，返回exp A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 0
if(n == 1) return poly(1, 1);
int t = (n + 1) >> 1;
poly res = exp(poly(A.begin(), A.begin() + t), t);
res = res * (poly(1, 1) - ln(res, n) + A);
res.resize(n); return res;
}
inline poly power(poly A, int k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
long long val = ksm(A[p], k), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(k * ln(A, n), n);
p = min(1ll * p * k, 1ll * n);
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly power(poly A, string k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
long long mod1 = 0, mod2 = 0;
bool zero = false; // 答案是否为 0
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
for(int i = 0; i < k.size(); i ++) {
if((mod1 * 10 + (k[i] - '0')) * p >= n) {zero = true; break; }
mod1 = (mod1 * 10 + (k[i] - '0')) % MOD;
mod2 = (mod2 * 10 + (k[i] - '0')) % (MOD - 1); // 指数对 \varphi(MOD) 取模
}
if(zero) return poly(n, 0);
long long val = ksm(A[p], mod2), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(mod1 * ln(A, n), n);
p = p * mod1;
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly sqrt1(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时请保证 A[0] = 1，可以稍加修改得到 A[0] 为完全平方数时的算法
int h = 1; while(h < n) h <<= 1;
A.resize(h, 0); poly res(1, 1);
for(int i = 2; i <= h; i <<= 1) {
res = (res * res + poly(A.begin(), A.begin() + i)) * inv(2 * res, i);
res.resize(i);
}
res.resize(n); return res;
}
// 二次剩余相关
long long w;
struct Complex {
int x, y;
Complex(int _x = 0, int _y = 0): x(_x), y(_y) {}
};
Complex operator * (const Complex &lhs, const Complex &rhs)
{return Complex(Plus(1ll * lhs.x * rhs.x % MOD, 1ll * lhs.y * rhs.y % MOD * w % MOD),
Plus(1ll * lhs.x * rhs.y % MOD, 1ll * lhs.y * rhs.x % MOD)); }
inline Complex complex_ksm(Complex A, int b) {
Complex r(1, 0);
for(; b; b >>= 1, A = A * A)
if(b & 1) r = r * A;
return r;
}
long long Cipolla(long long x) {
if (ksm(x,(MOD - 1) >> 1) == MOD - 1) return -1;
while(true) {
long long a = Rand() % MOD;
w = (a * a % MOD + MOD - x) % MOD;
if(ksm(w,(MOD - 1) >> 1) == MOD - 1) {
int res = complex_ksm(Complex(a, 1), (MOD + 1) >> 1).x;
return std::min(res, MOD - res); // 选用首项较小的解
}
}
}
inline poly sqrt(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时不需要保证 A[0] = 1
if(A.empty()) return A;
A.resize(n, 0);
long long val = A[0], mul = ksm(A[0], MOD - 2);
for(int i = 0; i < n; i ++) A[i] = A[i] * mul % MOD;
A = sqrt1(A, n);
val = Cipolla(val);
for(int i = 0; i < n; i ++) A[i] = A[i] * val % MOD;
return A;
}
inline void Rev(poly &A) {
// 反转系数
int n = (int)A.size() - 1;
for(int i = 0; 2 * i < n; i ++)
swap(A[i], A[n - i]);
}
inline pair<poly, poly> Div(poly A, poly B) {
// 多项式除法，返回的 std::pair 中第一维是商，第二维是余数
while(!A.empty() && A.back() == 0) A.pop_back();
while(!B.empty() && B.back() == 0) B.pop_back();
if(A.size() < B.size()) return {poly(1, 0), A};
int n = (int)A.size() - 1, m = (int)B.size() - 1;
Rev(A), Rev(B); poly p = A * inv(B, n - m + 1); p.resize(n - m + 1, 0);
Rev(A), Rev(B), Rev(p); poly r = A - B * p; r.resize(m, 0);
return {p, r};
}
inline void init() {make_iv(); }
} // namespace My_poly
using namespace My_poly;
int main() {
ios::sync_with_stdio(false), cin.tie(0);
init(); // poly 初始化
return 0;
}
遗憾的是，上面封装的多项式全家桶中，多项式exp的效率是普通方法的4倍，笔者尝试了优化但没能成功。
