多项式全家桶

发表于 更新于 本文字数： 11k 阅读时长 ≈ 10 分钟

默写课文。。。

导数公式：

原函数 导函数
$ y = C(C为常数) $ $ y’ = 0 $
$ y = a^x $ $ y’ = a^x \ln a $
$ y = e^x $ $ y’ = e^x $
$ y = x^n $ $ y’ = nx^{n - 1} $
$ y = \log_{a}x $ $ y’ = \frac{1}{x \ln a} $
$ y = \ln x $ $ y’ = \frac{1}{x} $
$ y = \sin x $ $ y’ = \cos x $
$ y = \cos x $ $ y’ = - \sin x $
$ y = u \pm v $ $ y’ = u’ \pm v’ $
$ y = uv $ $ y’ = u’v + uv’ $
$ y = \frac{u}{v} $ $ y’ = \frac{u’v - uv’}{v^2} $
$ y = f(g(x)) $ $ y’ = g’(x)f’(g(x)) $

完整代码：

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
inline int Max(const int &x, const int &y) { return x > y ? x : y; }
inline int Min(const int &x, const int &y) { return x < y ? x : y; }
inline int add(const int &x, const int &y) { return x + y < mod ? x + y : x + y - mod; }
inline int sub(const int &x, const int &y) { return x - y < 0 ? x - y + mod : x - y; }
inline int mul(const int &x, const int &y) { return (int)((ll)x * y % mod); }
inline int ksm(int x, int y = mod - 2) {
	int ss = 1; for(; y; y >>= 1, x = mul(x, x)) if(y & 1) ss = mul(ss, x); return ss;
}
namespace Poly {
	inline int Get(int x) { int ss = 1; for(; ss <= x; ss <<= 1); return ss; }
	void ntt(vector<int> &A, int lmt, int opt) {
		A.resize(lmt + 5);
		for(int i = 0, j = 0; i < lmt; i++) {
			if(i > j) swap(A[i], A[j]);
			for(int k = lmt >> 1; (j ^= k) < k; k >>= 1);
		}
		vector<int> w(lmt >> 1);
		for(int mid = 1; mid < lmt; mid <<= 1) {
			w[0] = 1;
			int w0 = ksm(opt == 1 ? 3 : (mod + 1) / 3, (mod - 1) / mid / 2);
			for(int j = 1; j < mid; j++) w[j] = mul(w[j - 1], w0);
			for(int R = mid << 1, j = 0; j < lmt; j += R)
				for(int k = 0; k < mid; k++) {
					int x = A[j + k], y = mul(w[k], A[j + mid + k]);
					A[j + k] = add(x, y), A[j + mid + k] = sub(x, y);
				}
		}
		if(opt == -1)
			for(int i = 0, inv = ksm(lmt); i < lmt; i++) A[i] = mul(A[i], inv);
	}
	vector<int> Add(const vector<int> &a, const vector<int> &b) {
		vector<int> res(Max(a.size(), b.size()));
		for(int i = 0; i < (int)res.size(); i++) {
			if(i < (int)a.size()) res[i] = add(res[i], a[i]);
			if(i < (int)b.size()) res[i] = add(res[i], b[i]);
		}
		return res;
	}
	vector<int> Mul(const vector<int> &a, const vector<int> &b) {
		vector<int> A = a, B = b; int lmt = Get(a.size() + b.size() - 2);
		ntt(A, lmt, 1), ntt(B, lmt, 1);
		for(int i = 0; i < lmt; i++) A[i] = mul(A[i], B[i]);
		ntt(A, lmt, -1); return A.resize(a.size() + b.size() - 1), A;
	}
	vector<int> Inv(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) return res.pb(ksm(A[0])), res;
		res = Inv(A, (sz + 1) / 2);
		vector<int> tmp(A.begin(), A.begin() + sz);
		int lmt = Get(sz * 2 - 2);
		ntt(tmp, lmt, 1), ntt(res, lmt, 1);
		for(int i = 0; i < lmt; i++) res[i] = mul(sub(2, mul(res[i], tmp[i])), res[i]);
		ntt(res, lmt, -1); return res.resize(sz), res;
	}
	void Div(const vector<int> &A, const vector<int> &B, vector<int> &D, vector<int> &R) {
		if(B.size() > A.size()) return (void)(D.clear(), D.pb(0), R = A);
		vector<int> a = A, b = B, iB; int n = A.size(), m = B.size();
		reverse(a.begin(), a.end()), reverse(b.begin(), b.end());
		b.resize(n - m + 1), iB = Inv(b, n - m + 1);
		D = Mul(a, iB), D.resize(n - m + 1), reverse(D.begin(), D.end());
		R = Mul(B, D);
		for(int i = 0; i < m - 1; i++) R[i] = (mod + A[i] - R[i]) % mod;
		R.resize(m - 1);
	}
	vector<int> Ln(const vector<int> &A) {
		vector<int> f, g = Inv(A); f.resize(A.size());
		for(int i = 1; i < (int)f.size(); i++) f[i - 1] = mul(i, A[i]);
		f[f.size() - 1] = 0, f = Mul(f, g), f.resize(A.size());
		for(int i = (int)f.size() - 1; i > 0; i--) f[i] = mul(f[i - 1], ksm(i));
		f[0] = 0; return f;
	}
	vector<int> Exp(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) return res.pb(1), res;
		res = Exp(A, (sz + 1) / 2), res.resize(sz); vector<int> tmp = Ln(res);
		for(int i = 0; i < (int)tmp.size(); i++) tmp[i] = add(sub(0, tmp[i]), A[i]);
		tmp[0] = add(1, tmp[0]), res = Mul(res, tmp);
		return res.resize(sz), res;
	}
	vector<int> Ksm(const vector<int> &A, const int &K) {
		int sz = A.size(), low = 0; vector<int> ta, res;
		for(int i = 0; i < sz; i++) if(A[i]) { low = i; break; }
		for(int i = low; i < sz; i++) ta.pb(A[i]);
		int mu = ta[0], inv = ksm(mu), Mu = ksm(mu, K);
		for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * inv % mod;
		ta = Ln(ta);
		for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * K % mod;
		ta = Exp(ta);
		for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * Mu % mod;
		int cnt = 0;
		for(int i = 0; i < low; i++)
			for(int j = 1; j <= K; j++) { ++cnt, res.pb(0); if(cnt > sz) break; }
		for(int i = 0; i < (int)ta.size(); i++) res.pb(ta[i]);
		return res.resize(sz), res;
	}
	vector<int> Ksm(const vector<int> &A, const vector<int> &K) {
		int mu = 0;
		for(int i = 0; i < (int)K.size(); i++) mu = add(mul(mu, 10), K[i]);
		return Ksm(A, mu);
	}
	vector<int> Sqrt(const vector<int> &A, int sz = -1) {
		if(sz == -1) sz = A.size();
		vector<int> res; if(sz == 1) { return res.pb(1), res; }
		res = Sqrt(A, (sz + 1) / 2); vector<int> tmp(A.begin(), A.begin() + sz);
		res.resize(sz), tmp = Mul(tmp, Inv(res));
		for(int i = 0; i < sz; i++) res[i] = 1ll * (res[i] + tmp[i]) * (mod + 1) / 2 % mod;
		return res;
	}
	vector<int> tt[4000010], res;
	vector<int> Evaluate_init(int u, int l, int r, const vector<int> &a) {
		vector<int> ttt;
		if(l >= r) { ttt.pb(mod - a[l]), ttt.pb(1); return tt[u] = ttt; }
		int mid = (l + r) >> 1;
		Evaluate_init(u * 2, l, mid, a), Evaluate_init(u * 2 + 1, mid + 1, r, a);
		return tt[u] = Mul(tt[u * 2], tt[u * 2 + 1]);
	}
	void Evaluate(int u, int l, int r, const vector<int> &f, const vector<int> &a) {
		if(r - l + 1 <= 512) {
			for(int k = l; k <= r; k++) {
				int w = 1, x = 0;
				for(int i = 0; i < (int)f.size(); i++) x = add(x, mul(w, f[i])), w = mul(w, a[k]);
				res.pb(x);
			}
			return ;
		}
		int mid = (l + r) >> 1; vector<int> tmp; Div(f, tt[u], tmp, tmp);
		Evaluate(u * 2, l, mid, tmp, a), Evaluate(u * 2 + 1, mid + 1, r, tmp, a);
	}
	vector<int> Evaluation(const vector<int> &f, const vector<int> &a) {
		Evaluate_init(1, 0, a.size() - 1, a);
		res.clear(), Evaluate(1, 0, a.size() - 1, f, a); return res;
	}
}

NTT（Mul&ntt）

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
inline int Get(int x) { int ss = 1; for(; ss <= x; ss <<= 1); return ss; }
void ntt(vector<int> &A, int lmt, int opt) {
	A.resize(lmt + 5);
	for(int i = 0, j = 0; i < lmt; i++) {
		if(i > j) swap(A[i], A[j]);
		for(int k = lmt >> 1; (j ^= k) < k; k >>= 1);
	}
	vector<int> w(lmt >> 1);
	for(int mid = 1; mid < lmt; mid <<= 1) {
		w[0] = 1;
		int w0 = ksm(opt == 1 ? 3 : (mod + 1) / 3, (mod - 1) / mid / 2);
		for(int j = 1; j < mid; j++) w[j] = mul(w[j - 1], w0);
		for(int R = mid << 1, j = 0; j < lmt; j += R)
			for(int k = 0; k < mid; k++) {
				int x = A[j + k], y = mul(w[k], A[j + mid + k]);
				A[j + k] = add(x, y), A[j + mid + k] = sub(x, y);
			}
	}
	if(opt == -1)
		for(int i = 0, inv = ksm(lmt); i < lmt; i++) A[i] = mul(A[i], inv);
}
vector<int> Mul(const vector<int> &a, const vector<int> &b) {
	vector<int> A = a, B = b; int lmt = Get(a.size() + b.size() - 2);
	ntt(A, lmt, 1), ntt(B, lmt, 1);
	for(int i = 0; i < lmt; i++) A[i] = mul(A[i], B[i]);
	ntt(A, lmt, -1); return A.resize(a.size() + b.size() - 1), A;
}

牛顿迭代：

$$
\begin{aligned}
\large{F(x) = F_0(x) - \frac{G(F_0(x))}{G’(F_0(x))}}\\
\large{G(F(x)) = 0}\\
\large{G’(F(x)) = \frac{dG}{dF}}\\
\end{aligned}
$$

推导：

$$
\begin{aligned}
y &= G’(F_0(x))(x - F_0(x)) + G(F_0(x))\\
0 &= G’(F_0(x))(F(x) - F_0(x)) + G(F_0(x))\\
F(x) &= F_0(x) - \frac{G(F_0(x))}{G’(F_0(x))}\\
\end{aligned}
$$

求逆（Inv）：

$$
\begin{aligned}
F \times G’ &= 1 (\mod x^{\lceil\frac{n}{2}\rceil})\\
F \times G &= 1 (\mod x^n)\\
F(G - G’) &= 0 (\mod x^{\lceil\frac{n}{2}\rceil})\\
G - G’ &= 0 (\mod x^{\lceil\frac{n}{2}\rceil})\\
(G - G’)^2 &= 0 (\mod x^n)\\
G^2 -2GG’ + G’^2 &= 0 (\mod x^n)\\
F(G^2 -2GG’ + G’^2) &= 0 (\mod x^n)\\
G &= 2G’ - FG’^2 (\mod x^n)\\
\end{aligned}
$$

1
2
3
4
5
6
7
8
9
10
11
vector<int> Inv(const vector<int> &A, int sz = -1) {
	if(sz == -1) sz = A.size();
	vector<int> res; if(sz == 1) return res.pb(ksm(A[0])), res;
	res = Inv(A, (sz + 1) / 2);
	vector<int> tmp(A.begin(), A.begin() + sz);
	int lmt = Get(sz * 2 - 2);
	ntt(tmp, lmt, 1), ntt(res, lmt, 1);
	for(int i = 0; i < lmt; i++)
		res[i] = mul(sub(2, mul(res[i], tmp[i])), res[i]);
	ntt(res, lmt, -1); return res.resize(sz), res;
}

求Ln（Ln）：

$$
\begin{aligned}
\ln’(F(x)) &= \frac{F’(x)}{F(x)}\\
\ln(F(x)) &= \int \ln’(F(x))\\
\end{aligned}
$$

1
2
3
4
5
6
7
vector<int> Ln(const vector<int> &A) {
	vector<int> f, g = Inv(A); f.resize(A.size());
	for(int i = 1; i < f.size(); i++) f[i - 1] = mul(i, A[i]);
	f[f.size() - 1] = 0, f = Mul(f, g), f.resize(A.size());
	for(int i = f.size() - 1; i > 0; i--) f[i] = mul(f[i - 1], ksm(i));
	f[0] = 0; return f;
}

求Exp（Exp）：

$$
\begin{aligned}
考虑牛顿迭代。\\
F(x) &= e^{a(x)}, \ln(F(x)) = a(x), G(F(x)) = ln(F(x)) - a(x) .\\
F(x) &= F_0(x) - \frac{G(F(x))}{G’(F(x))}\\
&= F_0(x) - \frac{\ln(F_0(x)) - a(x)}{\frac{1}{F_0(x)}}\\
&= F_0(x) - F_0(x)(\ln(F_0(x)) - a(x))\\
&= F_0(x)(1 - \ln(F_0(x)) + a(x))
\end{aligned}
$$

1
2
3
4
5
6
7
8
vector<int> Exp(const vector<int> &A, int sz = -1) {
	if(sz == -1) sz = A.size();
	vector<int> res; if(sz == 1) return res.pb(1), res;
	res = Exp(A, (sz + 1) / 2), res.resize(sz); vector<int> tmp = Ln(res);
	for(int i = 0; i < tmp.size(); i++) tmp[i] = add(sub(0, tmp[i]), A[i]);
	tmp[0] = add(1, tmp[0]), res = Mul(res, tmp);
	return res.resize(sz), res;
}

多项式快速幂（Ksm）：

$$
\begin{aligned}
(f(x))^{w} &= \exp(w\ln(f(x)))
\end{aligned}
$$

这种做法也可以用来多项式开根。。。 

1
2
3
4
5
6
vector<int> Ksm(const vector<int> &A, const vector<int> &K) {
	vector<int> res = Ln(A); int k = 0;
	for(int i = 0; i < K.size(); i++) k = add(mul(k, 10), K[i]);
	for(int i = 0; i < res.size(); i++) res[i]= mul(res[i], k);
	res = Exp(res), res.resize(A.size()); return res;
}

Update: 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
vector<int> Ksm(const vector<int> &A, const int &K) {
	int sz = A.size(), low = 0; vector<int> ta, res;
	for(int i = 0; i < sz; i++) if(A[i]) { low = i; break; }
	for(int i = low; i < sz; i++) ta.pb(A[i]);
	int mu = ta[0], inv = ksm(mu), Mu = ksm(mu, K);
	for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * inv % mod;
	ta = Ln(ta);
	for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * K % mod;
	ta = Exp(ta);
	for(int i = 0; i < (int)ta.size(); i++) ta[i] = 1ll * ta[i] * Mu % mod;
	int cnt = 0;
	for(int i = 0; i < low; i++)
		for(int j = 1; j <= K; j++) { ++cnt, res.pb(0); if(cnt > sz) break; }
	for(int i = 0; i < (int)ta.size(); i++) res.pb(ta[i]);
	return res.resize(sz), res;
}
vector<int> Ksm(const vector<int> &A, const vector<int> &K) {
	int mu = 0;
	for(int i = 0; i < (int)K.size(); i++) mu = add(mul(mu, 10), K[i]);
	return Ksm(A, mu);
}

开方（Sqrt）：

$$
\begin{aligned}
牛顿迭代。\\
G(F(x)) &= F^2(x) - a(x) = 0\\
F(x) &= F_0(x) - \frac{G(F_0(x))}{G’(F_0(x))}\\
&= F_0(x) - \frac{F_0^2(x) - a(x)}{2F_0(x)}\\
&= \frac{F_0^2(x) + a(x)}{2F_0(x)}\\
\end{aligned}
$$

1
2
3
4
5
6
7
8
vector<int> Sqrt(const vector<int> &A, int sz = -1) {
	if(sz == -1) sz = A.size();
	vector<int> res; if(sz == 1) { return res.pb(1), res; }
	res = Sqrt(A, (sz + 1) / 2); vector<int> tmp(A.begin(), A.begin() + sz);
	res.resize(sz), tmp = Mul(tmp, Inv(res));
	for(int i = 0; i < sz; i++) res[i] = 1ll * (res[i] + tmp[i]) * (mod + 1) / 2 % mod;
	return res;
}