本页面包含多项式常见的初等函数操作.具体而言,本页面包含如下内容:
- 多项式求逆
- 多项式开方
- 多项式除法
- 多项式取模
- 多项式指数函数
- 多项式对数函数
- 多项式三角函数
- 多项式反三角函数
初等函数与非初等函数
初等函数的定义如下1:
若域 F 中存在映射 u→∂u 满足:
- ∂(u+v)=∂u+∂v
- ∂(uv)=u∂v+v∂u
则称这个域为 微分域.
若微分域 F 上的函数 u 满足以下的任意一条条件,则称该函数 u 为初等函数:
- u 是 F 上的代数函数.
- u 是 F 上的指数性函数,即存在 a∈F 使得 ∂u=u∂a.
- u 是 F 上的对数性函数,即存在 a∈F 使得 ∂u=a∂a.
以下是常见的初等函数:
-
代数函数:存在有限次多项式 P 使得 P(f(x))=0 的函数 f(x),如 2x+1,x,(1+x2)−1,∣x∣.
-
指数函数
-
对数函数
-
三角函数
-
反三角函数
-
双曲函数
-
反双曲函数
-
以上函数的复合,如:
1+x2etanxsin(1+ln2x)
−iln(x+i1−x2)
以下是常见的非初等函数:
-
误差函数:
erf(x):=π2∫0xexp(−t2)dt
多项式求逆
给定多项式 f(x),求 f−1(x).
解法
倍增法
首先,易知
[x0]f−1(x)=([x0]f(x))−1
假设现在已经求出了 f(x) 在模 x⌈2n⌉ 意义下的逆元 f0−1(x).
有:
f(x)f0−1(x)f(x)f−1(x)f−1(x)−f0−1(x)≡1≡1≡0(modx⌈2n⌉)(modx⌈2n⌉)(modx⌈2n⌉)
两边平方可得:
f−2(x)−2f−1(x)f0−1(x)+f0−2(x)≡0(modxn)
两边同乘 f(x) 并移项可得:
f−1(x)≡f0−1(x)(2−f(x)f0−1(x))(modxn)
递归计算即可.
时间复杂度
T(n)=T(2n)+O(nlogn)=O(nlogn)
Newton's Method
参见 Newton's Method.
Graeffe 法
欲求 f−1(x)modx2n 考虑
f−1(x)modx2n=f(−x)(f(x)f(−x))−1modx2n=f(−x)g−1(x2)modx2n
只需求出 g−1(x)modxn 即可还原出 g−1(x2)modx2n 因为 f(x)f(−x) 是偶函数,时间复杂度同上.
代码
多项式求逆
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void polyinv(const poly &h, const int n, poly &f) {
static poly_t inv_t;
std::fill(f, f + n + n, 0);
f[0] = fpow(h[0], mod - 2);
for (int t = 2; t <= n; t <<= 1) {
const int t2 = t << 1;
std::copy(h, h + t, inv_t);
std::fill(inv_t + t, inv_t + t2, 0);
DFT(f, t2);
DFT(inv_t, t2);
for (int i = 0; i != t2; ++i)
f[i] = (i64)f[i] * sub(2, (i64)f[i] * inv_t[i] % mod) % mod;
IDFT(f, t2);
std::fill(f + t, f + t2, 0);
}
}
例题
- 有标号简单无向连通图计数:「POJ 1737」Connected Graph
多项式开方
给定多项式 g(x),求 f(x),满足:
f2(x)≡g(x)(modxn)
解法
倍增法
首先讨论 [x0]g(x) 不为 0 的情况.
易知:
[x0]f(x)=[x0]g(x)
若 [x0]g(x) 没有平方根,则多项式 g(x) 没有平方根.
[x0]g(x) 可能有多个平方根,选取不同的根会求出不同的 f(x).
假设现在已经求出了 g(x) 在模 x⌈2n⌉ 意义下的平方根 f0(x),则有:
f02(x)f02(x)−g(x)(f02(x)−g(x))2(f02(x)+g(x))2(2f0(x)f02(x)+g(x))22f0(x)f02(x)+g(x)2−1f0(x)+2−1f0−1(x)g(x)≡g(x)≡0≡0≡4f02(x)g(x)≡g(x)≡f(x)≡f(x)(modx⌈2n⌉)(modx⌈2n⌉)(modxn)(modxn)(modxn)(modxn)(modxn)
倍增计算即可.
时间复杂度
T(n)=T(2n)+O(nlogn)=O(nlogn)
还有一种常数较小的写法就是在倍增维护 f(x) 的时候同时维护 f−1(x) 而不是每次都求逆.
当 [x0]g(x)=1 时,可能需要使用二次剩余来计算 [x0]f(x).
上述方法需要知道 f0(x) 的逆,所以常数项不能为 0.
若 [x0]g(x)=0,则将 g(x) 分解成 xkh(x),其中 [x0]h(x)=0.
-
若 k 是奇数,则 g(x) 没有平方根.
-
若 k 是偶数,则求出 h(x) 的平方根 h(x),然后得到 f(x)≡xk/2h(x)(modxn).
洛谷模板题 P5205【模板】多项式开根 参考代码
#include <algorithm>
#include <iostream>
using namespace std;
constexpr int MAXN = 1 << 20, mod = 998244353;
int a[MAXN], b[MAXN], g[MAXN], gg[MAXN];
int qpow(int x, int y) {
int ans = 1;
while (y) {
if (y & 1) {
ans = (long long)1 * ans * x % mod;
}
x = (long long)1 * x * x % mod;
y >>= 1;
}
return ans;
}
int inv2 = qpow(2, mod - 2);
void change(int *f, int len) {
for (int i = 1, j = len / 2; i < len - 1; i++) {
if (i < j) {
swap(f[i], f[j]);
}
int k = len / 2;
while (j >= k) {
j -= k;
k /= 2;
}
if (j < k) {
j += k;
}
}
}
void NTT(int *f, int len, int type) {
change(f, len);
for (int q = 2; q <= len; q <<= 1) {
int nxt = qpow(3, (mod - 1) / q);
for (int i = 0; i < len; i += q) {
int w = 1;
for (int k = i; k < i + (q >> 1); k++) {
int x = f[k];
int y = (long long)1 * w * f[k + (q / 2)] % mod;
f[k] = (x + y) % mod;
f[k + (q / 2)] = (x - y + mod) % mod;
w = (long long)1 * w * nxt % mod;
}
}
}
if (type == -1) {
reverse(f + 1, f + len);
int iv = qpow(len, mod - 2);
for (int i = 0; i < len; i++) {
f[i] = (long long)1 * f[i] * iv % mod;
}
}
}
void inv(int deg, int *f, int *h) {
if (deg == 1) {
h[0] = qpow(f[0], mod - 2);
return;
}
inv((deg + 1) >> 1, f, h);
int len = 1;
while (len < deg * 2) {
len *= 2;
}
copy(f, f + deg, gg);
fill(gg + deg, gg + len, 0);
NTT(gg, len, 1);
NTT(h, len, 1);
for (int i = 0; i < len; i++) {
h[i] = (long long)1 * (2 - (long long)1 * gg[i] * h[i] % mod + mod) % mod *
h[i] % mod;
}
NTT(h, len, -1);
fill(h + deg, h + len, 0);
}
int n, t[MAXN];
void sqrt(int deg, int *f, int *h) {
if (deg == 1) {
h[0] = 1;
return;
}
sqrt((deg + 1) >> 1, f, h);
int len = 1;
while (len < deg * 2) {
len *= 2;
}
fill(g, g + len, 0);
inv(deg, h, g);
copy(f, f + deg, t);
fill(t + deg, t + len, 0);
NTT(t, len, 1);
NTT(g, len, 1);
NTT(h, len, 1);
for (int i = 0; i < len; i++) {
h[i] = (long long)1 * inv2 *
((long long)1 * h[i] % mod + (long long)1 * g[i] * t[i] % mod) % mod;
}
NTT(h, len, -1);
fill(h + deg, h + len, 0);
}
int main() {
cin >> n;
for (int i = 0; i < n; i++) {
cin >> a[i];
}
sqrt(n, a, b);
for (int i = 0; i < n; i++) {
cout << b[i] << ' ';
}
return 0;
}
Newton's Method
参见 Newton's Method.
例题
- 「Codeforces Round #250」E. The Child and Binary Tree
多项式除法 & 取模
给定多项式 f(x),g(x),求 g(x) 除 f(x) 的商 Q(x) 和余数 R(x).
解法
发现若能消除 R(x) 的影响则可直接 多项式求逆 解决.
考虑构造变换
fR(x)=xdegff(x1)
观察可知其实质为反转 f(x) 的系数.
设 n=degf,m=degg.
将 f(x)=Q(x)g(x)+R(x) 中的 x 替换成 x1 并将其两边都乘上 xn,得到:
xnf(x1)fR(x)=xn−mQ(x1)xmg(x1)+xn−m+1xm−1R(x1)=QR(x)gR(x)+xn−m+1RR(x)
注意到上式中 RR(x) 的系数为 xn−m+1,则将其放到模 xn−m+1 意义下即可消除 RR(x) 带来的影响.
又因 QR(x) 的次数为 (n−m)<(n−m+1),故 QR(x) 不会受到影响.
则:
fR(x)≡QR(x)gR(x)(modxn−m+1)
使用多项式求逆即可求出 Q(x),将其反代即可得到 R(x).
时间复杂度 O(nlogn).
多项式对数函数 & 指数函数
给定多项式 f(x),求模 xn 意义下的 lnf(x) 与 expf(x).
解法
普通方法
多项式对数函数
首先,对于多项式 f(x),若 lnf(x) 存在,则由其 定义,其必须满足:
[x0]f(x)=1
对 lnf(x) 求导再积分,可得:
dxdlnf(x)lnf(x)≡f(x)f′(x)≡∫dlnf(x)≡∫f(x)f′(x)dx(modxn)(modxn)
多项式的求导,积分时间复杂度为 O(n),求逆时间复杂度为 O(nlogn),故多项式求 ln 时间复杂度 O(nlogn).
多项式指数函数
首先,对于多项式 f(x),若 expf(x) 存在,则其必须满足:
[x0]f(x)=0
否则 expf(x) 的常数项不收敛.
对 expf(x) 求导,可得:
dxdexpf(x)≡expf(x)f′(x)(modxn)
比较两边系数可得:
[xn−1]dxdexpf(x)=i=0∑n−1([xi]expf(x))([xn−i−1]f′(x))
n[xn]expf(x)=i=0∑n−1([xi]expf(x))((n−i)[xn−i]f(x))
使用分治 FFT 即可解决.
时间复杂度 O(nlog2n).
Newton's Method
使用 Newton's Method 即可在 O(nlogn) 的时间复杂度内解决多项式 exp.
代码
多项式 ln/exp
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void derivative(const poly &h, const int n, poly &f) {
for (int i = 1; i != n; ++i) f[i - 1] = (i64)h[i] * i % mod;
f[n - 1] = 0;
}
void integrate(const poly &h, const int n, poly &f) {
for (int i = n - 1; i; --i) f[i] = (i64)h[i - 1] * inv[i] % mod;
f[0] = 0;
}
void polyln(const poly &h, const int n, poly &f) {
assert(h[0] == 1);
static poly_t ln_t;
const int t = n << 1;
derivative(h, n, ln_t);
std::fill(ln_t + n, ln_t + t, 0);
polyinv(h, n, f);
DFT(ln_t, t);
DFT(f, t);
for (int i = 0; i != t; ++i) ln_t[i] = (i64)ln_t[i] * f[i] % mod;
IDFT(ln_t, t);
integrate(ln_t, n, f);
}
void polyexp(const poly &h, const int n, poly &f) {
assert(h[0] == 0);
static poly_t exp_t;
std::fill(f, f + n + n, 0);
f[0] = 1;
for (int t = 2; t <= n; t <<= 1) {
const int t2 = t << 1;
polyln(f, t, exp_t);
exp_t[0] = sub(pls(h[0], 1), exp_t[0]);
for (int i = 1; i != t; ++i) exp_t[i] = sub(h[i], exp_t[i]);
std::fill(exp_t + t, exp_t + t2, 0);
DFT(f, t2);
DFT(exp_t, t2);
for (int i = 0; i != t2; ++i) f[i] = (i64)f[i] * exp_t[i] % mod;
IDFT(f, t2);
std::fill(f + t, f + t2, 0);
}
}
例题
-
计算 fk(x)
普通做法为多项式快速幂,时间复杂度 O(nlognlogk).
当 [x0]f(x)=1 时,有:
fk(x)=exp(klnf(x))
当 [x0]f(x)=1 时,设 f(x) 的最低次项为 fixi,则:
fk(x)=fikxikexp(klnfixif(x))
时间复杂度 O(nlogn).
多项式三角函数
给定多项式 f(x),求模 xn 意义下的 sinf(x),cosf(x) 与 tanf(x).
解法
首先由 Euler's formula (eix=cosx+isinx) 可以得到 三角函数的另一个表达式:
sinxcosx=2ieix−e−ix=2eix+e−ix
那么代入 f(x) 就有:
sinf(x)cosf(x)=2iexp(if(x))−exp(−if(x))=2exp(if(x))+exp(−if(x))
直接按上述表达式编写程序即可得到模 xn 意义下的 sinf(x) 与 cosf(x).再由 tanf(x)=cosf(x)sinf(x) 可求得 tanf(x).
代码
多项式三角函数
注意到我们是在 Z998244353 上做 NTT,那么相应地,虚数单位 i 应该被换成 86583718 或 911660635:
⟹i=−1≡998244352(mod998244353)ori≡86583718(mod998244353)ori≡911660635(mod998244353)
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
constexpr int imgunit = 86583718;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void polytri(const poly &h, const int n, poly &sin_t, poly &cos_t) {
assert(h[0] == 0);
static poly_t tri1_t, tri2_t;
for (int i = 0; i != n; ++i) tri2_t[i] = (i64)h[i] * imgunit % mod;
polyexp(tri2_t, n, tri1_t);
polyinv(tri1_t, n, tri2_t);
if (sin_t != nullptr) {
const int invi = fpow(pls(imgunit, imgunit), mod - 2);
for (int i = 0; i != n; ++i)
sin_t[i] = (i64)(tri1_t[i] - tri2_t[i] + mod) * invi % mod;
}
if (cos_t != nullptr) {
for (int i = 0; i != n; ++i) cos_t[i] = div2(pls(tri1_t[i], tri2_t[i]));
}
}
多项式反三角函数
给定多项式 f(x),求模 xn 意义下的 arcsinf(x),arccosf(x) 与 arctanf(x).
解法
仿照求多项式 ln 的方法,对反三角函数求导再积分可得:
dxdarcsinxarcsinxdxdarccosxarccosxdxdarctanxarctanx=1−x21=∫1−x21dx=−1−x21=−∫1−x21dx=1+x21=∫1+x21dx
那么代入 f(x) 就有:
dxdarcsinf(x)arcsinf(x)dxdarccosf(x)arccosf(x)dxdarctanf(x)arctanf(x)=1−f2(x)f′(x)=∫1−f2(x)f′(x)dx=−1−f2(x)f′(x)=−∫1−f2(x)f′(x)dx=1+f2(x)f′(x)=∫1+f2(x)f′(x)dx
直接按式子求就可以了.
代码
多项式反三角函数
constexpr int MAXN = 262144;
constexpr int mod = 998244353;
using i64 = long long;
using poly_t = int[MAXN];
using poly = int *const;
void derivative(const poly &h, const int n, poly &f) {
for (int i = 1; i != n; ++i) f[i - 1] = (i64)h[i] * i % mod;
f[n - 1] = 0;
}
void integrate(const poly &h, const int n, poly &f) {
for (int i = n - 1; i; --i) f[i] = (i64)h[i - 1] * inv[i] % mod;
f[0] = 0;
}
void polyarcsin(const poly &h, const int n, poly &f) {
static poly_t arcsin_t;
const int t = n << 1;
std::copy(h, h + n, arcsin_t);
std::fill(arcsin_t + n, arcsin_t + t, 0);
DFT(arcsin_t, t);
for (int i = 0; i != t; ++i) arcsin_t[i] = sqr(arcsin_t[i]);
IDFT(arcsin_t, t);
arcsin_t[0] = sub(1, arcsin_t[0]);
for (int i = 1; i != n; ++i)
arcsin_t[i] = arcsin_t[i] ? mod - arcsin_t[i] : 0;
polysqrt(arcsin_t, n, f);
polyinv(f, n, arcsin_t);
derivative(h, n, f);
DFT(f, t);
DFT(arcsin_t, t);
for (int i = 0; i != t; ++i) arcsin_t[i] = (i64)f[i] * arcsin_t[i] % mod;
IDFT(arcsin_t, t);
integrate(arcsin_t, n, f);
}
void polyarccos(const poly &h, const int n, poly &f) {
polyarcsin(h, n, f);
for (int i = 0; i != n; ++i) f[i] = f[i] ? mod - f[i] : 0;
}
void polyarctan(const poly &h, const int n, poly &f) {
static poly_t arctan_t;
const int t = n << 1;
std::copy(h, h + n, arctan_t);
std::fill(arctan_t + n, arctan_t + t, 0);
DFT(arctan_t, t);
for (int i = 0; i != t; ++i) arctan_t[i] = sqr(arctan_t[i]);
IDFT(arctan_t, t);
inc(arctan_t[0], 1);
std::fill(arctan_t + n, arctan_t + t, 0);
polyinv(arctan_t, n, f);
derivative(h, n, arctan_t);
DFT(f, t);
DFT(arctan_t, t);
for (int i = 0; i != t; ++i) arctan_t[i] = (i64)f[i] * arctan_t[i] % mod;
IDFT(arctan_t, t);
integrate(arctan_t, n, f);
}
参考资料与链接