从零开始的多项式笔记(四 - 多项式基本操作)

本文是基于 oi-wiki 整理的,大部分是照抄,修改了少量错误;和本博客、oi-wiki 原文一样遵循 CC BY-SA 4.0 协议

[TOC]


多项式乘法逆

problem

给定多项式 $f(x)$,求 $f^{-1}(x)$。

solution - 倍增法

首先,易知

假设现在已经求出了 $f(x)$ 在模 $x^{ \lfloor n / 2 \rfloor}$ 意义下的逆元 $f_0^{-1}(x)$,有

递归计算即可。复杂度 $O(n \log n)$。

code - 递归实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void Inv(int N, int *a, int *b) {
if (N == 1) {
b[0] = Pow(a[0], mod - 2, mod);
return;
}
// 清空
Inv((N + 1) >> 1, a, b);

int lim, l, *c = T;
for (lim = 1, l = -1; lim <= N << 1; lim <<= 1, ++l);
for (int i = 0; i < lim; ++i)
to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);

memcpy(c, a, sizeof(int) * N);
memset(c + N, 0, sizeof(int) * (lim - N));
NTT(lim, b, 1), NTT(lim, c, 1);
for (int i = 0; i < lim; ++i)
b[i] = (2ll - 1ll * b[i] * c[i] % mod + mod) % mod * b[i] % mod;
NTT(lim, b, -1);
memset(b + N, 0, sizeof(int) * (lim - N));
}

code - 迭代实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void Inv(int N, int *a, int *b) {
int *c = T1;
memset(b, 0, sizeof(int) * (N << 1));
b[0] = Pow(a[0], mod - 2, mod);
for (int m = 2, lim = 4, l = 1; m <= N << 1; m <<= 1, lim <<= 1, ++l) {
for (int i = 0; i < lim; ++i)
to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);

memcpy(c, a, sizeof(int) * m);
memset(c + m, 0, sizeof(int) * (lim - m));
NTT(lim, b, 1), NTT(lim, c, 1);
for (int i = 0; i < lim; ++i)
b[i] = (2ll - 1ll * b[i] * c[i] % mod + mod) % mod * b[i] % mod;
NTT(lim, b, -1);
memset(b + m, 0, sizeof(int) * (lim - m));
}
}


多项式开方

problem

给定多项式 $f(x)$,求 $g(x)$ 满足

solution - 倍增法

假设现在已经求出了 $f(x)$ 在模 $x^{ \lfloor n / 2 \rfloor}$ 意义下的平方根 $g_0(x)$,有

递归计算即可。复杂度 $O(n \log n)$。

在不保证保证 $[x_0] f(x) = 1$ 时需要整二次剩余算 $[x_0] g(x)$。

code - 迭代实现

巨大巨大的常数,跑得是论文鸽十倍多慢,我爬罢

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
void Sqrt(int N, int *a, int *b) {
int *c = T1, *d = T1;
memset(b, 0, sizeof(int) * (N << 1));
b[0] = 1;
for (int m = 2, lim = 4, l = 1; m < N << 1; m <<= 1, lim <<= 1, ++l) {
memset(c, 0, sizeof(int) * lim);
Inv(m, b, c);

// for (int i = 0; i < lim; ++i)
// to[i] = (to[i >> 1] >> 1) | ((i & 1) << l);
// 这两行可以不要是因为在求逆过程的最后也会计算 to[]

memcpy(d, a, sizeof(int) * m);
memset(d + m, 0, sizeof(int) * (lim - m));
NTT(lim, b, 1), NTT(lim, c, 1), NTT(lim, d, 1);
for (int i = 0; i < lim; ++i)
b[i] = (b[i] + 1ll * c[i] * d[i]) % mod * inv2 % mod;
NTT(lim, b, 0);
memset(b + m, 0, sizeof(int) * (lim - m));
}
}

多项式除法

problem

给定度数为 $n$ 的多项式 $f(x)$ 和度数为 $m \ (m \le n)$ 的多项式 $g(x)$,求 $Q(x), R(x)$ 满足

且 $\operatorname{deg} Q \le n - m, \operatorname{deg} R < m$。

solution

发现比较麻烦的是 $R(x)$,消掉 $R(x)$ 后就是一个多项式求逆。考虑构造变换

其实质为翻转 $F(x)$ 的系数。

将 $f(x) = Q(x) g(x) + R(x)$ 中的 $x$ 替换成 $1 / x$ 并将等式两边同乘 $x^n$ 得

上式中 $R^R(x)$ 的系数为 $x^{n - m + 1}$,则将其放到模 $x^{n - m + 1}$ 意义下即可消除 $R^R(x)$ 的影响;又 $\operatorname{deg} Q^R = n - m < n - m + 1$,所以 $Q^R(x)$ 不会受影响。则有

多项式求逆求出 $Q(x)$,反代得到 $R(x)$。总复杂度 $O(n \log n)$。

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
void Div(int N, int M, int *f, int *g, int *Q, int *R) {
int lim, m = N - M + 1, *t = T1, *s = T2;

for (lim = 1; lim < m << 1; lim <<= 1);
std::reverse_copy(g, g + M, t); // t <- Rev(g)
Inv(m, t, s); // s <- Inv(Rev(g))
memset(s + m, 0, sizeof(int) * (lim - m)); // s <- Rev(g) mod x^(deg f - deg g + 1)
std::reverse_copy(f, f + N, t); // t <- Rev(f)
memset(t + m, 0, sizeof(int) * (lim - m)); // t <- Rev(f) mod x^(deg f - deg g + 1)

NTT(lim, t, 1), NTT(lim, s, 1);
for (int i = 0; i < lim; ++i)
t[i] = 1ll * t[i] * s[i] % mod; // t <- Rev(Q)
NTT(lim, t, 0);

for (lim = 1; lim < N; lim <<= 1);
memset(t + m, 0, sizeof(int) * (lim - m));
std::reverse(t, t + m); // t <- Q
memcpy(Q, t, sizeof(int) * m);
memcpy(s, g, sizeof(int) * M); // s <- g
memset(s + M, 0, sizeof(int) * (lim - M));

NTT(lim, t, 1), NTT(lim, s, 1);
for (int i = 0; i < lim; ++i)
t[i] = 1ll * t[i] * s[i] % mod; // t <- Q * s
NTT(lim, t, 0);
for (int i = 0; i < M; ++i)
R[i] = (f[i] - t[i] + mod) % mod;
}


多项式求导 & 积分

1
2
3
4
5
6
7
8
9
10
11
void Deri(int N, int *f, int *g) { // [0, N) -> [0, N - 1)
for (int i = 1; i < N; ++i)
g[i - 1] = 1ll * i * f[i] % mod;
g[N - 1] = 0;
}

void Int(int N, int *f, int *g) { // [0, N) -> [0, N + 1)
for (int i = N; i; --i)
g[i] = 1ll * f[i - 1] * inv[i] % mod;
g[0] = 0;
}

多项式对数函数

problem

给定多项式 $f(x)$,求 $\ln(f(x))$。保证 $[x^0] f(x) = 1$。

solution

对于多项式 $f(x)$,若 $\ln f(x)$ 存在,其必须满足

对 $\ln f(x)$ 求导再积分可得

多项式求导、积分复杂度为 $O(n)$,求逆和乘法复杂度为 $O(n \log n)$,总复杂度 $O(n \log n)$。

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void Ln(int N, int *f, int *g) {
int lim, *t = T2;

Inv(N, f, g); // g <- Inv(f)
Deri(N, f, t); // t <- f'
memset(t + N, 0, sizeof(int) * N);

for (lim = 1; lim < N << 1; lim <<= 1);
NTT(lim, g, 1), NTT(lim, t, 1);
for (int i = 0; i < lim; ++i)
g[i] = 1ll * g[i] * t[i] % mod; // g <- f' * Inv(f)
NTT(lim, g, 0);

Int(N, g, g);
}


多项式指数函数

problem

给定多项式 $f(x)$,求 $\exp(f(x))$。保证 $[x^0] f(x) = 0$。

solution

对于多项式 $f(x)$,若 $\ln f(x)$ 存在,其必须满足

否则 $\exp f(x)$ 的常数项不收敛。

但是直接做这玩意这玩意要分治 FFT 我不会啊,直接去牛顿迭代里学不用分治 FFT 复杂度还低的算法吧:

复杂度 $O(n \log n)$。

code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void Exp(int N, int *f, int *g) {
int *t = T3;
g[0] = 1;
for (int m = 2, lim = 4; m < N << 1; m <<= 1, lim <<= 1) {
Ln(m, g, t); // t <- Ln(f_0)
for (int i = 0; i < m; ++i)
t[i] = (f[i] - t[i] + mod) % mod;
t[0] = (t[0] + 1) % mod; // t <- 1 - Ln(f_0) + h

NTT(lim, g, 1), NTT(lim, t, 1);
for (int i = 0; i < lim; ++i)
g[i] = 1ll * g[i] * t[i] % mod;
NTT(lim, g, 0);

memset(g + m, 0, sizeof(int) * (lim - m));
}
}


多项式幂函数

problem

给定多项式 $f(x)$ 和整数 $k$,求 $f^k(x)$。

solution

显然可以直接套快速幂,复杂度 $O(n \log n \log k)$。有时候 $k$ 会非常大,需要一个更优秀的算法。

设 $f(x)$ 的最低次项为 $f_i x^i$,则

特别地,当 $[x^0] f(x) = 1$ 时,有

复杂度 $O(n \log n)$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
void Pow(int N, int K, int *f, int *g) { // 保证 f[0] = 1
int *t = T4;

Ln(N, f, t);
for (int i = 0; i < N; ++i)
t[i] = 1ll * t[i] * K % mod;
Exp(N, t, g);
}

void Pow(int N, int K, int K_, int *f, int *g) {
// K 为指数模 mod 的值,K_ 为指数模 phi(mod) 的值
// 下文的 sp 代表指数是否大于等于 mod
int low, *t = T4;

for (low = 0; !f[low]; ++low);
if (1ll * low * K >= N || (low && sp)) return;

int len_ = low * K, len = N - len_;
int tk = ::Pow(f[low], K_), invT = ::Pow(f[low]);
for (int i = low; i < N; ++i)
t[i] = 1ll * f[i] * invT % mod;
Ln(len, t + low, g);
for (int i = 0; i < len; ++i)
g[i] = 1ll * g[i] * K % mod;
Exp(len, g, t);
memset(g, 0, sizeof(int) * len_);
for (int i = 0; i < len; ++i)
g[i + len_] = 1ll * t[i] * tk % mod;
}
文章作者: fa_555
文章链接: https://blog.fa555.tech/2020/Polynomial-4/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 fa_555's Blog