mydcwfy's Blog

多项式汇总

Posted on
Symbols count in article: 18k Reading time ≈ 40 mins.

大杂烩。建议先学 FFT 和 NTT。

也可以来这里找代码。

注意清空数组！

会不断更新。

这里汇总了所有我所学过的多项式的板子。

既是存代码，也是简单的讲解。

首先甩一个 NTT 在这里。

前置知识：FFT。

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
typedef long long ll;

ll qpow(ll a, ll k)
{
    ll res = 1;
    while (k)
    {
        if (k & 1) res = res * a % Mod;
        a = a * a % Mod;
        k >>= 1;
    }
    return res;
}

void poly_rev(int bit)
{
    for (int i = 1; i < (1 << bit); ++ i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
}

void NTT(ll *a, int bit, int inv)
{
    poly_rev(bit);
    int tot = 1 << bit;
    for (int i = 0; i < tot; ++ i)
        if (rev[i] < i) swap(a[rev[i]], a[i]);
    for (int mid = 1; mid < tot; mid <<= 1)
    {
        ll ak = qpow(3, (Mod - 1) / (mid << 1));
        if (inv == -1) ak = qpow(ak, Mod - 2);
        for (int i = 0; i < tot; i += (mid << 1))
        {
            ll now = 1;
            for (int j = 0; j < mid; ++ j, now = now * ak % Mod)
            {
                ll x = a[i + j], y = a[i + j + mid] * now % Mod;
                a[i + j] = (x + y) % Mod, a[i + j + mid] = (x - y + Mod) % Mod;
            }
        }
    }
    if (inv == 1) return;
    ll Inv = qpow(tot, Mod - 2);
    for (int i = 0; i < tot; ++ i) a[i] = a[i] * Inv % Mod;
}

1. 多项式求逆

可以考虑倍增。

默认 $\dfrac{n}{2}=\lfloor\dfrac{n+1}{2}\rfloor$

假设我们已经求出了 $h(x)\equiv f(x)^{-1}\pmod{x^{\frac{n}{2}}}$。

我们需要求 $g(x)\equiv f(x)^{-1}\pmod{x^n}$。

很明显，可以得到：$g(x)-h(x)\equiv 0\pmod {x^{\frac{n}{2} } }$。

平方一下，就是 $g^2(x)-2g(x)h(x)+h^2(x)\equiv 0\pmod {x^n}$。

有一个 $g^2(x)$，我们考虑搞掉一个，所以乘一个 $f(x)$，也就是 $g(x)-2h(x)+h^2(x)f(x)\equiv 0\pmod {x^n}$。

于是：$g(x)\equiv h(x)(2-h(x)f(x))\pmod {x^n}$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void poly_inv(ll *a, ll *b, int len)
{
    if (len == 1)
    {
        b[0] = qpow(a[0], Mod - 2);
        return;
    }
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    poly_inv(a, b, (len + 1) >> 1);
    static ll c[N];
    for (int i = (len + 1) >> 1; i < tot; ++ i) b[i] = 0;
    for (int i = 0; i < tot; ++ i)
        if (i < len) c[i] = a[i];
        else c[i] = 0;
    NTT(c, bit, 1), NTT(b, bit, 1);
    for (int i = 0; i < tot; ++ i)
        b[i] = (2 - b[i] * c[i] % Mod + Mod) % Mod * b[i] % Mod;
    NTT(b, bit, -1);
    for (int i = len; i < tot; ++ i) b[i] = 0;
}

2. 多项式 ln

考虑：$g(x)=\ln f(x)$。

同时求导：$g’(x)=\ln f(x)$。

看做复合函数，$h(x)=\ln x$，那么 $\ln f(x) = h(f(x)) = \dfrac{1}{f(x)} \cdot f’(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
void poly_dao(ll *a, ll *b, int len)
{
    for (int i = 1; i < len; ++ i) b[i - 1] = a[i] * i % Mod;
    b[len - 1] = 0;
}

void poly_ji(ll *a, ll *b, int len)
{
    for (int i = 1; i < len; ++ i) b[i] = a[i - 1] * qpow(i, Mod - 2) % Mod;
    b[0] = 0;
}

void poly_ln(ll *a, ll *b, int len)
{
    static ll c[N], d[N];
    poly_dao(a, c, len), poly_inv(a, d, len);
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = len; i < tot; ++ i) c[i] = d[i] = 0;
    NTT(c, bit, 1), NTT(d, bit, 1);
    for (int i = 0; i < tot; ++ i) c[i] = c[i] * d[i] % Mod;
    NTT(c, bit, -1);
    for (int i = len; i < tot; ++ i) c[i] = 0;
    poly_ji(c, b, len);
    for (int i = len; i < tot; ++ i) b[i] = 0;
}

3. 多项式 exp

这个需要牛顿迭代。

假设我们要找到一个函数的零点，我们当前找到的值为 $x_0$，怎样才能逼近更优解呢？

画一个图。

看到上面的函数图像，我们发现，可以直接对该点进行求导，然后交到 $x$ 轴的位置，会更靠近答案。

于是我们可以得到式子：

精度就可以不断提高。

回到本题，经过我不会的证明，可以得到，单次牛顿迭代后，答案可以由 $\pmod{x ^ { \frac{n} {2} } }$ 变为 $\pmod{ x^n }$。

那么，假设我们已经确定了 $h(x) = e ^ {f(x)} \pmod {x ^ {\frac {n}{2} } }$，要求 $g(x) = e ^{f(x)} \pmod {x ^ n}$。

我们考虑构造 $A(x) = \ln g(x) - f(x) = 0$，求 $g(x)$，现在我们已经找到了一个近似的 $h(x)$。

带入上面的式子，就是：

直接 NTT 就可以了。注意保证 $f(0) = 0$，也就是常数项为 1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void poly_exp(ll *a, ll *b, int len)
{
    if (len == 1)
    {
        b[0] = 1;
        return;
    }
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    poly_exp(a, b, (len + 1) >> 1);
    static ll d[N];
    for (int i = 0; i < tot; ++ i) d[i] = 0;
    poly_ln(b, d, len);
    for (int i = 0; i < len; ++ i) d[i] = (a[i] - d[i] + Mod) % Mod;
    for (int i = len; i < tot; ++ i) d[i] = 0;
    d[0] ++;
    NTT(d, bit, 1), NTT(b, bit, 1);
    for (int i = 0; i < tot; ++ i)
        b[i] = b[i] * d[i] % Mod;
    NTT(b, bit, -1);
    for (int i = len; i < tot; ++ i) b[i] = 0;
}

4. 多项式快速幂

比较简单。

唯一需要注意的是如果 $[x^0]A(x) \not = 1$ 的时候，而 $\ln$ 的要求就是 $[x^0]A(x) = 1$，我们需要将最低非 0 项移到最低点。所以要注意移位的问题和全部是 0 的情况。

这是保证 $[x^0]A(x) = 1$ 的情况的写法。

1
2
3
4
5
6
7
void poly_qpow(LL *a, LL *b, int len, int k)
{
    static LL c[N];
    poly_ln(a, c, len);
    for (int i = 0; i < len; ++ i) c[i] = c[i] * k % Mod;
    poly_exp(c, b, len);
}

这是 Luogu 模板题毒瘤写法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void poly_qpow(const ll *a, ll *b, int len, ll k1, ll k2)//k1 = k % Mod, k2 = k % (Mod - 1), flag = (k >= Mod)
{
    int t = 0;
    while (!a[t] && t < len) t ++;
    if (1ll * t * k1 >= len || (t && flag)) return;
    static ll c[N], d[N];
    ll Inv = qpow(a[t], Mod - 2), z = a[t];
    len -= t;
    for (int i = 0; i < len; ++ i) d[i] = a[i + t] * Inv % Mod;
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = 0; i < tot; ++ i) c[i] = 0;
    for (int i = len; i < tot; ++ i) d[i] = 0;
    poly_ln(d, c, len);
    for (int i = 0; i < len; ++ i) c[i] = c[i] * k1 % Mod;
    for (int i = len; i < tot; ++ i) c[i] = 0;
    poly_exp(c, b, len);
    len += t;
    for (int i = len - 1; i >= t * k1; -- i) b[i] = b[i - t * k1] * qpow(z, k2) % Mod;
    for (int i = t * k1 - 1; ~i; -- i) b[i] = 0;
}

5. 多项式开根

是一个简单的推导：

一般 $[x^0]A(x) \not= 0$，所以我们直接找到 $\sqrt {[x^0]A(x)}$，然后 $pow$ 之前每一项都除以 $[x^0]A(x)$，然后 $pow$，最后将每一项都乘上 $\sqrt{[x^0]A(x)}$ 就是了。至于如何在模意义下求根号，我们使用 Cipolla 算法。

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
namespace Cipolla{
    struct Complex{
        ll x, y;
    };
    int w;

    Complex mul(Complex a, Complex b)
    {
        return {(a.x * b.x % Mod + a.y * b.y % Mod * w % Mod) % Mod, (a.x * b.y % Mod + b.x * a.y % Mod) % Mod};
    }

    Complex qpow_img(Complex a, int k)
    {
        Complex res = {1, 0};
        while (k)
        {
            if (k & 1) res = mul(res, a);
            a = mul(a, a);
            k >>= 1;
        }
        return res;
    }

    ll get_sqrt(ll x)
    {
        ll a;
        if (qpow(x, (Mod - 1) >> 1) == Mod - 1) return -1;
        while (1)
        {
            a = rand() % Mod;
            w = (a * a % Mod - x + Mod) % Mod;
            if (qpow(w, (Mod - 1) >> 1) == Mod - 1) break;
        }
        ll ans = qpow_img({a, 1}, (Mod + 1) >> 1).x;
        return min(ans, Mod - ans);
    }
};
using Cipolla::get_sqrt;

void poly_sqrt(ll *f, ll *g, int len)
{
    ll sq = get_sqrt(f[0]), inv = qpow(f[0], Mod - 2);
    for (int i = 0; i < len; ++ i) f[i] = f[i] * inv % Mod;
    poly_qpow(f, g, len, (Mod + 1) >> 1, (Mod + 1) >> 1);
    for (int i = 0; i < len; ++ i) g[i] = g[i] * sq % Mod;
}

6. 多项式除法

我们假设这样一个函数：

表示将 $F(x)$ 的系数翻转过来。可以发现 $F(x) = x^n F(\dfrac 1x)$。

然后，我们对这个式子进行推导（假设是 $F(x) = D(x) * G(x) + R(x)$）：

然后，如果我们将这一个看做一个 $\bmod x^{n - m}$ 的话，那么 $R_R(x) * x ^{n - m}$ 就没有用了。

然后，$D_R(x)$ 又是 $n - m$ 次的，所以我们在 $\bmod x ^{n - m}$ 下用求逆求出 $D_R(x)$ 就可以了。至于要求出 $R(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
void poly_div(LL *a, LL *b, LL *c, LL *r, int n, int m)
{
    int del = n - m + 1;
    static LL d[N], e[N];
    int bit = 0;
    while ((1 << bit) < (n << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = 0; i < tot; ++ i)
        if (i < m) e[i] = b[m - i - 1];
        else e[i] = 0;
    poly_inv(e, d, del);
    for (int i = 0; i < tot; ++ i)
        if (i < n) e[i] = a[n - i - 1];
        else e[i] = 0;
    NTT(d, bit, 1), NTT(e, bit, 1);
    for (int i = 0; i < tot; ++ i)
        d[i] = d[i] * e[i] % Mod;
    NTT(d, bit, -1);
    for (int i = 0; i < del; ++ i) c[i] = d[del - i - 1];

    for (int i = 0; i < tot; ++ i)
        if (i < del) d[i] = c[i];
        else d[i] = 0;
    for (int i = 0; i < tot; ++ i)
        if (i < m) e[i] = b[i];
        else e[i] = 0;
    NTT(d, bit, 1), NTT(e, bit, 1);
    for (int i = 0; i < tot; ++ i) d[i] = d[i] * e[i] % Mod;
    NTT(d, bit, -1);
    for (int i = n; i < tot; ++ i) d[i] = 0;
    for (int i = 0; i < n; ++ i)
        r[i] = (a[i] - d[i] + Mod) % Mod;
}

7. 多项式三角函数

一个大名鼎鼎的公式：

带入 $-x$ 得：

（不会有人连三角函数的诱导公式都不知道吧，$\cos(-x) = \cos x, \sin (-x) = -\sin x$）。

于是我们直接解出来：

至于 $i$ 在模 998244353 下的取值，我们有 $i^4 = 1 = g^{p - 1}$，所以 $i = g^{\frac{p - 1}{4}}$，在 998244353 下是 86583718。

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
const int img = 86583718;

void poly_sin(ll *a, ll *b, int len)
{
    static ll c[N], d[N], e[N];
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = 0; i < len; ++ i) c[i] = a[i] * img % Mod;
    for (int i = len; i < tot; ++ i) c[i] = 0;
    poly_exp(c, d, len);
    for (int i = 0; i < len; ++ i) c[i] = Mod - a[i] * img % Mod;
    for (int i = len; i < tot; ++ i) c[i] = 0;
    poly_exp(c, e, len);
    ll invimg = qpow(2ll * img % Mod, Mod - 2);
    for (int i = 0; i < len; ++ i) b[i] = (d[i] - e[i] + Mod) % Mod * invimg % Mod;
}

void poly_cos(ll *a, ll *b, int len)
{
    static ll c[N], d[N], e[N];
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = 0; i < len; ++ i) c[i] = a[i] * img % Mod;
    for (int i = len; i < tot; ++ i) c[i] = 0;
    poly_exp(c, d, len);
    for (int i = 0; i < len; ++ i) c[i] = Mod - a[i] * img % Mod;
    poly_exp(c, e, len);
    ll inv = (Mod + 1) >> 1;
    for (int i = 0; i < len; ++ i) b[i] = (d[i] + e[i]) % Mod * inv % Mod;
}

8. 多项式反三角函数

先扔一个公式（注意 $\text{asin}$ 和 $\arcsin$ 似乎是一个东西）：

简单的证明一下（应该高数上有吧）：

其中第四步是复合函数的求导：$y = g(t), t = f(x)\Rightarrow y’ = g’(t) * f’(x)$。$\arctan$ 的略去了。

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
void poly_asin(LL *a, LL *b, int len)
{
    static LL c[N], d[N], e[N];
    poly_dao(a, c, len);
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = 0; i < tot; ++ i)
        if (i < len) d[i] = a[i];
        else d[i] = 0;
    NTT(d, bit, 1);
    for (int i = 0; i < tot; ++ i) d[i] = d[i] * d[i] % Mod;
    NTT(d, bit, -1);
    for (int i = 0; i < tot; ++ i)
        if (i < len) d[i] = (!i - d[i] + Mod) % Mod;
        else d[i] = 0;
    poly_sqrt(d, e, len);
    for (int i = 0; i < tot; ++ i) d[i] = 0;
    poly_inv(e, d, len);
    NTT(c, bit, 1), NTT(d, bit, 1);
    for (int i = 0; i < tot; ++ i) c[i] = c[i] * d[i] % Mod;
    NTT(c, bit, -1);
    poly_ji(c, b, len);
}

void poly_atan(LL *a, LL *b, int len)
{
    static LL c[N], d[N], e[N];
    poly_dao(a, c, len);
    int bit = 0;
    while ((1 << bit) < (len << 1)) bit ++;
    int tot = 1 << bit;
    for (int i = 0; i < tot; ++ i)
        if (i < len) d[i] = a[i];
        else d[i] = 0;
    NTT(d, bit, 1);
    for (int i = 0; i < tot; ++ i) d[i] = d[i] * d[i] % Mod;
    NTT(d, bit, -1);
    for (int i = 0; i < tot; ++ i)
        if (i < len) d[i] = !i + d[i];
        else d[i] = 0;
    poly_inv(d, e, len);
    NTT(c, bit, 1), NTT(e, bit, 1);
    for (int i = 0; i < tot; ++ i) c[i] = c[i] * e[i] % Mod;
    NTT(c, bit, -1);
    poly_ji(c, b, len);
}
点击有惊喜 
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
/// updated on 2022-03-18: used adj(x) to calculate mod, changed N's value
/// It is written by mydcwfy on 2021-12-30.
/// This is the template of Polynomial, including NTT, Polyinverse, PolyDerivate, PolyIntegrate, PolyLn, PolyExp, PolyQpow, PolySqrt, Polysin..., Polyasin...
/// The caculation is under the mod of 998244353.
/// You can use the namespace under your version.
namespace Polynomial {

const int N = 1 << 21 | 1;
const int Mod = 998244353, Img = 86583718;

int rev[N], Inv[N];
Cipolla cip(998244353);


/// Calculate the mod of x
/// please make sure that x is in [-Mod, Mod)
inline int adj(int x) { return x += x >> 31 & Mod; }

/// This is the background fuction of the Polynomial.
/// It caculates the k-pow of a.
LL qpow(LL a, LL k)
{
	LL res = 1;
	while (k)
	{
		if (k & 1) res = res * a % Mod;
		a = a * a % Mod;
		k >>= 1;
	}
	return res;
}

/// This is the background fuction of the Polynomial.
/// It caculates the reverse of NTT.
void Poly_rev(int bit)
{
	for (int i = 1; i < (1 << bit); ++ i)
		rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
}

/// This is the background fuction of the Polynomial.
/// It caculates the mininum number bit qualifying (2 ^ bit >= len * 2)
int Poly_bit(int len)
{
	int bit = 0;
	while ((1 << bit) < (len << 1)) bit ++;
	return bit;
}

/// This is the background fuction of the Polynomial.
/// It caculates the inverse of numbers ranging from [1, N).
void init()
{
	Inv[1] = 1;
	for (int i = 2; i < N; ++ i) Inv[i] = (Mod - Mod / i) * (LL)Inv[Mod % i] % Mod;
}
/// This is the background fuction of the Polynomial.
void NTT(LL *a, int bit, int inv)
{
	if (!Inv[1]) init();
	Poly_rev(bit);
	int tot = 1 << bit;
	for (int i = 0; i < tot; ++ i)
		if (i < rev[i]) swap(a[rev[i]], a[i]);
	for (int mid = 1, bit = 0; mid < tot; mid <<= 1, bit ++) {
		LL ak = qpow(inv == 1 ? 3 : (Mod + 1) / 3, (Mod - 1) >> (bit + 1));
		for (int i = 0; i < tot; i += mid << 1)
			for (int j = 0, now = 1; j < mid; ++ j, now = now * ak % Mod)
			{
				int x = a[i + j], y = a[i + j + mid] * now % Mod;
				a[i + j] = adj(x + y - Mod), a[i + j + mid] = adj(x - y);
			}
	}
	if (inv == 1) return;
	for (int i = 0; i < tot; ++ i) a[i] = a[i] * Inv[tot] % Mod;
}

/// Caculate the Polyinverse of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Inv(LL *a, LL *b, int len)
{
	if (len == 1)
	{
		b[0] = qpow(a[0], Mod - 2);
		return;
	}
	Poly_Inv(a, b, (len + 1) >> 1);
	static LL c[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = (len + 1) >> 1; i < tot; ++ i) b[i] = 0;
	for (int i = 0; i < tot; ++ i)
		if (i < len) c[i] = a[i];
		else c[i] = 0;
	NTT(c, bit, 1), NTT(b, bit, 1);
	for (int i = 0; i < tot; ++ i)
		b[i] = (2 - b[i] * c[i] % Mod + Mod) % Mod * b[i] % Mod;
	NTT(b, bit, -1);
	for (int i = len; i < tot; ++ i) b[i] = 0;
}

/// Caculate the PolyDerivate of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Dev(LL *a, LL *b, int len)
{
	for (int i = 1; i < len; ++ i) b[i - 1] = a[i] * i % Mod;
	b[len - 1] = 0;
}

/// Caculate the PolyIntegrate of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Int(LL *a, LL *b, int len)
{
	if (!Inv[1]) init();
	for (int i = 1; i < len; ++ i) b[i] = a[i - 1] * Inv[i] % Mod;
	b[0] = 0;
}

/// Caculate the PolyLn of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Ln(LL *a, LL *b, int len)
{
	static LL c[N], d[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = 0;
	Poly_Dev(a, c, len), Poly_Inv(a, d, len);
	NTT(c, bit, 1), NTT(d, bit, 1);
	for (int i = 0; i < tot; ++ i) c[i] = c[i] * d[i] % Mod;
	NTT(c, bit, -1);
	Poly_Int(c, b, len);
	for (int i = len; i < tot; ++ i) b[i] = 0;
}

/// Caculate the PolyExp of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Exp(LL *a, LL *b, int len)
{
	if (len == 1)
	{
		b[0] = 1;
		return;
	}
	Poly_Exp(a, b, (len + 1) >> 1);
	int bit = Poly_bit(len), tot = 1 << bit;
	static LL c[N];
	for (int i = 0; i < tot; ++ i) c[i] = 0;
	Poly_Ln(b, c, len);
	for (int i = 0; i < len; ++ i) c[i] = (!i + a[i] - c[i] + Mod) % Mod;
	NTT(b, bit, 1), NTT(c, bit, 1);
	for (int i = 0; i < tot; ++ i) b[i] = b[i] * c[i] % Mod;
	NTT(b, bit, -1);
	for (int i = len; i < tot; ++ i) b[i] = 0;
}

/// Caculate the PolyQpow of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
/// k is the pow number.
void Poly_Qpow(LL *a, LL *b, int len, LL k)
{
	static LL c[N];
	LL x = a[0], z = qpow(x, Mod - 2), mul = qpow(x, k % (Mod - 1));
	for (int i = 0; i < len; ++ i) a[i] = a[i] * z % Mod;
	for (int i = 0; i < len; ++ i) c[i] = 0;
	Poly_Ln(a, c, len);
	for (int i = 0; i < len; ++ i) c[i] = c[i] * (k % Mod) % Mod;
	Poly_Exp(c, b, len);
	for (int i = 0; i < len; ++ i) b[i] = b[i] * mul % Mod;
}

/// Caculate the PolySqrt of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Sqrt(LL *a, LL *b, int len)
{
	LL x = a[0], z = cip.Mod_Sqrt(x), cu = qpow(x, Mod - 2);
	for (int i = 0; i < len; ++ i) a[i] = a[i] * cu % Mod;
	Poly_Qpow(a, b, len, (Mod + 1) >> 1);
	for (int i = 0; i < len; ++ i) a[i] = a[i] * z % Mod;
}

/// Caculate the PolySin of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Sin(LL *a, LL *b, int len)
{
	static LL c[N], d[N], e[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = e[i] = 0;
	for (int i = 0; i < tot; ++ i)
		if (i < len) c[i] = a[i] * Img % Mod;
		else c[i] = 0;
	Poly_Exp(c, d, len);
	for (int i = 0; i < tot; ++ i)
		if (i < len) c[i] = Mod - a[i] * Img % Mod;
		else c[i] = 0;
	Poly_Exp(c, e, len);
	LL InvImg = qpow(2LL * Img % Mod, Mod - 2);
	for (int i = 0; i < len; ++ i) b[i] = (d[i] - e[i] + Mod) % Mod * InvImg % Mod;
}

/// Caculate the PolyCos of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Cos(LL *a, LL *b, int len)
{
	static LL c[N], d[N], e[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = e[i] = 0;
	for (int i = 0; i < len; ++ i) c[i] = a[i] * Img % Mod;
	Poly_Exp(c, d, len);
	for (int i = 0; i < tot; ++ i)
		if (i < len) c[i] = Mod - a[i] * Img % Mod;
		else c[i] = 0;
	Poly_Exp(c, e, len);
	for (int i = 0; i < len; ++ i) b[i] = (d[i] + e[i]) % Mod * Inv[2] % Mod;
}

/// Caculate the PolyTan of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Tan(LL *a, LL *b, int len)
{
	static LL c[N], d[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = 0;
	Poly_Sin(a, c, len), Poly_Cos(a, d, len), Poly_Inv(d, b, len);
	NTT(c, bit, 1), NTT(b, bit, 1);
	for (int i = 0; i < tot; ++ i) b[i] = c[i] * b[i] % Mod;
	NTT(b, bit, -1);
	for (int i = len; i < tot; ++ i) b[i] = 0;
}

/// Caculate the PolyCot of a.
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Cot(LL *a, LL *b, int len)
{
	static LL c[N], d[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = 0;
	Poly_Cos(a, c, len), Poly_Sin(a, d, len), Poly_Inv(d, b, len);
	NTT(c, bit, 1), NTT(b, bit, 1);
	for (int i = 0; i < tot; ++ i) b[i] = c[i] * b[i] % Mod;
	NTT(b, bit, -1);
	for (int i = len; i < tot; ++ i) b[i] = 0;
}

/// Caculate the PolyAsin of a(Arcsin).
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Asin(LL *a, LL *b, int len)
{
	static LL c[N], d[N], e[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = e[i] = 0;
	Poly_Dev(a, c, len);
	for (int i = 0; i < tot; ++ i)
		if (i < len) d[i] = a[i];
		else d[i] = 0;
	NTT(d, bit, 1);
	for (int i = 0; i < tot; ++ i) d[i] = d[i] * d[i] % Mod;
	NTT(d, bit, -1);
	for (int i = 0; i < tot; ++ i)
		if (i < len) d[i] = (!i - d[i] + Mod) % Mod;
		else d[i] = 0;
	Poly_Sqrt(d, e, len);
	for (int i = 0; i < tot; ++ i) d[i] = 0;
	Poly_Inv(e, d, len);
	NTT(c, bit, 1), NTT(d, bit, 1);
	for (int i = 0; i < tot; ++ i) c[i] = c[i] * d[i] % Mod;
	NTT(c, bit, -1);
	Poly_Int(c, b, len);
}

/// Caculate the PolyAtan of a(Arctan).
/// a is the original Polynomial.
/// b is the Polynomial to be caculated into.
/// len is the length of the Polynomial.
void Poly_Atan(LL *a, LL *b, int len)
{
	static LL c[N], d[N], e[N];
	int bit = Poly_bit(len), tot = 1 << bit;
	for (int i = 0; i < tot; ++ i) c[i] = d[i] = e[i] = 0;
	for (int i = 0; i < len; ++ i) c[i] = a[i];
	NTT(c, bit, 1);
	for (int i = 0; i < tot; ++ i) c[i] = c[i] * c[i] % Mod;
	NTT(c, bit, -1);
	for (int i = len; i < tot; ++ i) c[i] = 0;
	c[0] ++;
	Poly_Dev(a, d, len), Poly_Inv(c, e, len);
	NTT(d, bit, 1), NTT(e, bit, 1);
	for (int i = 0; i < tot; ++ i) d[i] = d[i] * e[i] % Mod;
	NTT(d, bit, -1);
	Poly_Int(d, b, len);
}

} // namespace Polynomial