多项式操作

多项式运算从入门到入坟

目录

• 多项式运算从入门到入坟
  • 1. 前言和前置知识
  • 2. 微积分基础知识
    • 2. 1 积分
    • 2. 2 导数
    • 2. 3 不定积分和微积分基本定理
    • 2. 4 基本函数求导公式
    • 2. 5 求导运算法则
    • 2. 6 常用积分公式
    • 2. 7 泰勒展开和泰勒级数
  • 3. 牛顿迭代
  • 4. 多项式求逆
  • 5. 多项式求导
  • 6. 多项式积分
  • 7. 多项式对数函数
  • 8. 多项式指数函数
  • 9. 多项式快速幂
  • 10. 多项式开方
  • 11. 多项式带余除法
  • 12. 多项式三角函数
  • 13. 多项式多点求值
  • 14. 多项式快速插值
  • 15. 拉格朗日反演
  • 16. 分治 NTT
  • 17. 完结撒花

1. 前言和前置知识

蛤？不会 FFT 你敢来这儿？

多项式是啥？详见 《普通初中课程标准实验教科书 数学 七年级 上册》（人民教育出版社）。

生成函数就是将数列 ({a_n}) 变成了这样一个多项式

[G(x) = sum_{k = 0}^{infty}a_kx^k ]

但是不要被那个 (infty) 吓到了，我们并不关心这个函数值是多少，我们只关心各项的系数。

而且我们在题里只需要知道前(n)项的系数就行了，本质上是用一个多项式来表示数列。

为什么要多项式表示呢？因为多项式有各种各样神仙操作，这就是我们下面要讲的各种东西。

为了充分利用多项式的神奇性质，我们需要知道这些东西：

• 高一基础数学知识
• 微积分基础
• NTT

本文并不会涉及到这些多项式操作的组合意义是什么，主要讲的就是板子怎么写，也就是最基础的部分。

但是其中有些毒瘤操作可能并不会用上，一般来讲掌握 (4 - 10) 节的内容便可以应对大部分的多项式题。

2. 微积分基础知识

别害怕，没什么难的。

主要讲一讲基础的部分，微积分是一门博大精深水深的学科，详细的部分大家会在大学学习，这里只讲一些皮毛。

2. 1 积分

记不记得高一时最开始学的物理？

没错就是运动学。

位移(x)，速度(v)，加速度(a)是什么关系？

我们来回顾一下匀加速实现运动：虽然速度 (v = at) ，但是我们怎么算位移呢？

由于速度是不断变化，我们并没有什么显然的式子来算位移。

但如果我们将速度看做一段段的匀速运动，我们就可以用式子 (x = vt) 来计算位移。

于是我们将速度分成 (Delta t) 长的一段段，总时间为 (t) ，我们就有了 (N = frac {t}{Delta t}) 个时间段。而对于第 (k) 个时间段，我们将这段时间内的速度近似的看成 (v_k = akDelta t) ，我们就可以将最后的位移式子近似的算成

[egin {align} x=&sum_{k = 1}^{N} v_kDelta t \ =& sum_{k = 1}^{N} akDelta t^2 \ =& aDelta t^2sum_{k = 1}^Nk \ =& aDelta t^2 frac {N(N+1)}{2} \ =& frac 1 2 aDelta t^2 left(left(frac {t}{Delta t}{} ight)^2 + left(frac {t}{Delta t}{} ight) ight) \ =& frac 1 2 at^2 + frac 1 2 at Delta t end {align} ]

然后当 $Delta t $ 趋于(0)的时候，$x =frac 1 2 at^2 $ 。怎么样，有没有很熟悉？还记得吗？

我们来思考一下， $Delta t $ 趋于(0)是啥意思？我们将 (t) 分成了 $ frac {t}{Delta t}$ 个时间段， $Delta t $ 趋于 (0) ，也就是说我们的段数趋近于无穷大。虽然 (Delta t) 等于 (0) 的时候我们会浮点数例外，但是我们却不能否认当 (Delta t) 越来越接近 (0) 的时候，这个式子的值永远不会小于 $frac 1 2 at^2 $ 。而这一道永远迈不过去的，真实存在的 “坎”，就是极限。这就是微积分中最基础的思想——极限。

以上过程实际上就是积分的过程，用更正规的表达方式来表示，就是：

[x = int_0^tvmathrm d t = frac 1 2 at^2 ]

在这里， (x) 和 (v) 都是一个关于时间 (t) 的函数。而 (d) 并没有实际含义，只是一个算子，表示我们的函数的自变量是 (t) 。

而积分的感性理解，就是求函数图像与 (x) 轴所围成的面积。

2. 2 导数

假如我们已经知道了位移的表达式，我们怎么知道速度和加速度的表达式呢？

我们来看速度的定义：(v = frac {Delta x} {Delta t}) ，也就是位移函数的斜率。当位移是一个匀加速运动的时候，我们可以这样容易地计算出速度。但如果位移的表达式是 (x =frac 1 2 at^2) 怎么办？

我们依旧将位移函数看成一段段的一次函数。我们取两个时间点 (t,t_1) ，满足 (t_1 = t + Delta t) 。于是我们的速度就变成了

[egin {align} v =& frac {Delta x} {Delta t} \ =& frac {x_{t_1}-x_t} {t_1-t} \ =& frac 1 2 a frac {t_1^2-t^2}{t_1 - t} \ =& frac 1 2 a frac {(t_1-t)(t_1+t)}{t_1 - t} \ =& frac 1 2 a (t_1 + t) end {align} ]

同样的，我们另 (Delta t) 趋于 (0) ，这样 (v = at) 。

当 (Delta t) 趋于 (0) 的时候，我们好像又浮点数例外了，然而还是那个极限的思想：虽然我们不能除 (0) ，但是这个极限却确实存在。所以这一点的斜率，我们就定义成这个极限。

（其实微积分中好多不合情理的东西最终都是个定义）。

于是我们也得到了导数的感性理解：函数图像上某一点的斜率。

2. 3 不定积分和微积分基本定理

刚才的积分是有上下界的，叫做定积分，它是一个值。

而如果没有上下界，叫做不定积分，它是一个函数。

不定积分的过程，是已知一个函数 (F(x)) 的导数 (f(x)) ，来求 (F(x )) 。我们称 (F(x)) 为 (f(x)) 的原函数。

（注意导数也是一个函数）。

我们接下来所说的多项式积分，是求不定积分，也就是求另一个多项式。

顺带提一下微积分基本定理：

[int_a^bf(x)mathrm dx = F(b) - F(a) = F(x)Large|_a^b ]

2. 4 基本函数求导公式

时间关系不再一一证明。

感兴趣的同学可以课下自己学习《微积分》还有b站3Blue1Brown的视频。

以下除了 (x) 之外的字母都是常数。 (y) 的导数表示成 (y') 。

[egin {align} c' &= 0 \ x^n {'} &= nx^{n-1} \ a^x {'} &= a^xln a \ e^x {'} &= e^x \ log_ax' &= frac {1} {x ln a} \ ln x ' &= frac {1} {x} \ sin x' &= cos x \ cos x' &= -sin x end {align} ]

2. 5 求导运算法则

(f(x),g(x)) 都表示关于 (x) 的函数。

[egin {align} (f(x) pm g(x))' &= f'(x) pm g'(x) \ (f(x)g(x))' &= f'(x)g(x) + f(x)g'(x) \ left(frac{f(x)}{g(x)} ight)' &= frac {f'(x)g(x) - f(x)g'(x)} {g(x)^2} \ (f(g(x)))' &= f'(g(x))g'(x) end {align} ]

2. 6 常用积分公式

只需要知道两个。

[egin {align} int x^nmathrm dx &= frac {1} {n+1}x^{n+1} + C (n e -1)\ int frac 1 x mathrm dx &= ln |x| + C end {align} ]

因为求导的话会消掉常数项，所以不定积分是要带一个常数 (C) 的，但是我们一般当成 (0)。

积分运算法则和 $ sum $ 一样（吧大概）。

2. 7 泰勒展开和泰勒级数

具体形式参见在美妙的数学王国中畅游。

已知 (f(x_0)) ，我们可以这样近似 (f( x )) （别问我为什么）

（ (f^{(n)}(x)) 表示 (f(x)) 的 (n) 阶导数，就是求导 (n) 次。）

[f(x) = sum_{n=0}^{infty}frac {f^{(n)}(x_0)(x-x_0)^n} {n!} ]

特别的，当 (x_0 = 0) 时，这个级数又叫麦克劳林级数。

3. 牛顿迭代

以下所有东西都能用这个东西推。

假如我们现在有一个多项式 (f(x)) 和一个方程 (g(f(x)) = 0) ，我们已知 (f(x)) 的前 (n) 项 (f_0(x)) ，要求 (f(x)) 的前 (2n) 项。

牛顿迭代就是用来解这个方程的。

我们有

[f(x) equiv f_0(x) mod x^n ]

我们对 (g(f(x))) 进行泰勒展开，得到

[egin {align} 0 &= g(f(x)) \ &= g(f_0(x)) + g'(f_0(x))(f(x)-f_0(x))+frac 1 2g''(f_0(x))(f(x)-f_0(x))^2+dots \ &equiv g(f_0(x)) + g'(f_0(x))(f(x)-f_0(x)) mod x^{2n} end {align} ]

因为 (f(x)-f_0(x)) 最低项的次数是 (n) ，所以泰勒展开中平方及以上的项都被模掉了。

化一化式子，我们就得到了

[f(x) equiv f_0(x) - frac{g(f_0(x))}{g'(f_0(x))} mod x^{2n} ]

每一次迭代，项数翻倍，复杂度

[T(n) = T(n/2) + O(n log n) = O(n log n) ]

4. 多项式求逆

以下模数默认 (998244353)。

已知 (A(x)) ，求 (B(x)) 满足 (A(x)B(x) equiv 1) 。

由于不明原因，这里用牛顿迭代解释不太科学。考虑这样构造。

[egin {align} A(x)B(x) &equiv 1 mod x^n\ A(x)B(x) - 1 &equiv 0 mod x^{n} \ (A(x)B(x) - 1) ^ 2 &equiv 0 mod x^{2n} \ A(x)(2B(x) - A(x)B(x)^2) &equiv 1 mod x^{2n} \ end {align} ]

5. 多项式求导

按照定义即可。

[A'(x) = sum_{k geq 1}ka_kx^{k-1} ]

6. 多项式积分

按照定义即可。

[int A(x) = sum_{k geq 0}frac {a_k x^{k+1}} {k+1} ]

7. 多项式对数函数

这个用不着迭代。

求 $B(x) = ln(A(x)) $ ，则 (B'(x)=frac {A'(x)}{A(x)})

然后多项式求导求逆加积分就好了。

8. 多项式指数函数

设 (B(x) = e^{A(x)})，

构造 (g(B(x)) = ln(B(x)) - A(x) = 0)。

直接套公式 (B(x) = B_0(x) + B_0(x)(A(x) - ln(B_0(x)))。

注意这里要保证 (A(x)) 的常数项 为 (0) 。不然你从哪儿找到一个 (e^c) 呢？

9. 多项式快速幂

先整体除常数项，算完后乘常数项的(k)次幂就好啦！

为什么要除呢？因为 (ln) 要保证常数项为 (1) 。

$A(x)^k = (e^{ln A(x)})^k = e^{kln A(x)} $ 。

10. 多项式开方

你可以选择牛顿迭代推式子，你也可以选择用上面的快速幂，也就是 ( ext{Inv}2) 次幂。

那么 ( ext{Inv2}) 是模多少意义下的呢？

其实是模 (MOD) 意义下而不是 (varphi(MOD) = MOD - 1) 意义下的。

虽然在指数上，可这只是一个记号而已。你并没有找到一个 (e) ，然后算 (e) 的多项式次幂。

你实际上算的是这样一个函数

[exp(F(x)) = sum_{k = 0}^infty frac{F(x)^i}{i!} ]

所以 (k) 并没有在指数上，自然也就不是模 (MOD - 1) 而是模 (MOD) 。

同理，如果你能处理常数项的问题，也就是解 (k) 次剩余，你就能给多项式开任意方根。

11. 多项式带余除法

这和求逆不太一样，要求给出一个商多项式和一个余多项式。

大概的思想就是将余多项式搞没，然后算出商，最后再算出来余多项式。

这里设

[egin{align} F(x) = sum_{i = 0}^nf_ix^i \ G(x) = sum_{i = 0}^mg_ix^i \ end{align} ]

然后要求你算出 (F(x) = G(x) *Q(x) + P(x)) 。

这里假设 (n geq m) , (Q(x)) 的次数为 (n - m) ，(P(x)) 的次数为 (m - 1) 。

考虑这样一个东西

[F_R(x) = x^nF(frac{1}{x}) ]

我们发现 (F_R(x)) 实际上就是把 (F(x)) 的系数翻转了。

那我们尝试推一下式子。

[egin {align} F(x) &= Q(x) * G(x) +P(x) \ F(frac 1 x) &= Q(frac 1 x) * G(frac 1 x) + P(frac 1 x) \ x^nF(frac 1 x) &= x^{n-m}Q(frac 1 x) * x^mG(frac 1 x) + x^{n-m+1}x^{m-1}P(frac 1 x) \ x^nF(frac 1 x) &= x^{n-m}Q(frac 1 x) * x^mG(frac 1 x) + x^{n-m+1}x^{m-1}P(frac 1 x) mod x^{n-m+1} \ x^nF(frac 1 x) &= x^{n-m}Q(frac 1 x) * x^mG(frac 1 x) mod x^{n-m+1} \ F_R(x) &= Q_R(x) * G_R(x) mod x^{n - m + 1} \ end {align} ]

然后我们普通的求逆求出来 (Q_R(x)) 之后，std::reverse 就能得到 (Q(x)) 。

而

[P(x) = F(x) - Q(x)*G(x) ]

我们也就可以得到 (P(x)) 。

12. 多项式三角函数

(sin) 和 (cos) 。

这个东西真的一点用都没有。但是我们就是能算。

首先根据泰勒展开我们有

[egin{align} exp(ix) = cos(x) + isin(x) \ exp(-ix) = cos(x) - isin(x) end{align} ]

（其实这就是欧拉公式的由来）。

所以

[egin{align} sin(x) = frac{exp(ix) - exp(-ix)}{2i} \ cos(x) = frac{exp(ix) + exp(-ix)}{2} end{align} ]

给读进来的多项式乘个 (i) ，(exp) 一下加加减减就出来了。

那么哪里来的 (i) 呢？

如果你像我一样傻，你会暴力试所有数找一个平方等于 (998244352) 的玩意。

然而 (i) 不就是个 (4) 次单位根吗？现成的啊。

qpow(3, (MOD - 1) / 4)

没了。

(arcsin) 和 (arctan)。

这玩意比上面那俩还没用。不过更好写。

根据高数内容，

[arcsin(x)' = frac{1}{sqrt{1-x^2}} \ arctan(x)' = frac{1}{1+x^2} \ ]

然后

[egin{align} arcsin(A(x)) = int frac{A(x)'}{sqrt{1-A(x)^2}} \ arctan(A(x)) = int frac{A(x)'}{1+A(x)^2} \ end {align} ]

没了。

13. 多项式多点求值

多点求值，就是给你一个多项式 (A(x)) 和一堆 (x) 坐标让你求对应的 (y) 。

构造多项式

[F_{l,r}(x) = prod_{i=l}^r(x - x_i) ]

然后让 (A(x)) 模 (F_{l,r}) ，得到一个余式 (R(x)) ，即

[A(x) = Q(x) F_{l,r}(x) + R(x) ]

然后发现当 (x) 取 (x_l) 到 (x_r) 中的值时，原多项式的值和那个余式的值一样。于是只用处理那个余式。

然而暴力模每个 (F_{l,l}) 复杂度并不对。于是我们可以先模 (F_{l,r}) ，然后把模剩下的多项式分别模 (F_{l,mid},F_{mid+1,r}) ，递归处理。每次取模时，余式的次数会变为模多项式的次数减一，于是每次次数减半，复杂度 (O(nlog^2n)) ，模到叶子上时正好就剩一个常数了，也就是那个对应的点值。

实现时分治加 NTT 求 (F) ，然后用一个线段树一样的数组结构存下来每一个 (F_{l,r}) 。

Poly mem[MAXN];
inline void Divide(int u, int l, int r, const Poly& A) {
  if (l == r) {
    mem[u] = Poly({MOD - A[l], 1});
    return;
  }
  int mid = (l + r) >> 1;
  Divide(u << 1, l, mid, A), Divide(u << 1 | 1, mid + 1, r, A);
  mem[u] = mem[u << 1] * mem[u << 1 | 1];
}

inline void Eval_recur(int u, int l, int r, const Poly& A, Poly& res) {
  if (l == r) return res.push_back(A[0]);
  int mid = (l + r) >> 1;
  Eval_recur(u << 1, l, mid, Div(A, mem[u << 1]).second, res);
  Eval_recur(u << 1 | 1, mid + 1, r, Div(A, mem[u << 1 | 1]).second, res);
}

inline Poly Multieval(const Poly& A, const Poly& B) {
  Divide(1, 0, B.size() - 1, B);
  Poly res;
  Eval_recur(1, 0, B.size() - 1, Div(A, mem[1]).second, res);
  return res;
}

14. 多项式快速插值

快速插值，就是 NTT 优化拉格朗日插值。

先看拉格朗日插值的式子

[F(x)=sum_{i=1}^n{frac{prod_{{j} eq{i}}{({x}-{x_j})}}{prod_{{j} eq{i}}{({x_i}-{x_j})}}{y_i}} ]

分母是个常数，单独和 (y) 提出来。

[F(x)=sum_{i=1}^n{frac{y_i}{prod_{{j} eq{i}}{({x_i}-{x_j})}}prod_{{j} eq{i}}{({x}-{x_j})}} ]

如果我们手头有一个

[G(x) = prod_{}{({x}-{x_i})} ]

那么我们就需要求出来这个多项式

[frac{G(x)}{(x - x_i)} ]

在 (x_i) 处的取值。

肯定不能每次都除啊，复杂度原地爆炸。

我们考虑直接把 (x_i) 代进去，就成了 (frac{0}{0}) 。

这时直接用洛必达法则求极限就是正确值。

什么是洛必达法则？

当你需要求 (frac{F(x)}{G(x)}) 在某一点 (x_0) 的极限的时候，很不巧，变成了 (frac 0 0) 或者 (frac infty infty) 这样的形式。

我们就给分子分母同时求导，一直导到分子分母不是 (frac infty infty) 和 (frac 0 0) 为止。这个时候极限就是 (frac{F(x_0)}{G(x_0)}) 。

然后我们给上面那个式子洛必达一下，得到 (G'(x_i)) 。

直接上多点求值就能把所有的 (frac{y_i}{prod_{{j} eq{i}}{({x_i}-{x_j})}}) 求出来。

后面的那个 (prod_{{j} eq{i}}{({x}-{x_j})}) 分治求就行了。

具体可以看代码实现。用到了一些上面多点求值的函数。

inline Poly Insv_recur(int u, int l, int r, const Poly& A) {
  if (l == r) return Poly({A[l]});
  int mid = (l + r) >> 1;
  Poly a = Insv_recur(u << 1, l, mid, A) * mem[u << 1 | 1];
  Poly b = Insv_recur(u << 1 | 1, mid + 1, r, A) * mem[u << 1];
  return a + b;
}

inline Poly Fastinsv(const Poly& X, const Poly& Y) {
  Divide(1, 0, X.size() - 1, X);
  Poly A = Derivate(mem[1]), B;
  Eval_recur(1, 0, X.size() - 1, A, B);
  for (size_t i = 0; i < X.size(); ++i)
    B[i] = mul(qpow(B[i], MOD - 2), Y[i]);
  return Insv_recur(1, 0, X.size() - 1, B);
}

讲完了。

15. 拉格朗日反演

又是这拉格朗日。

这里我们需要引入一个 “复合逆” 的概念。

“复合逆” 和反函数差不多。当 (F(G(x)) = 1) 的时候， (F(x)) 和 (G(x)) 互为复合逆。

我们现在的任务就是求一个多项式的复合逆。

然而不幸的是，我们并没有复杂度低于 (O(n^2)) 的求复合逆算法。

但是我们可以快速求得某一项的取值。

限于作者水平不能给出证明。这里只给出计算方法。

[[x^n] g(x) = frac 1 n [x^{n-1}]left(frac {x}{f(x)} ight)^n ]

而且我们还有广义拉格朗日反演

[[x^n]h(g(x))=frac 1 n[x^{n-1}]h'(x)left(frac x {f(x)} ight)^n ]

16. 分治 NTT

当卷积变成一种“我卷我自己”的形式的时候，我们可以考虑使用 CDQ 分治来求。

就是转移式成了这个样子

[f_i = sum_{j = 1}^i f_{i - j}g_j ]

我们就可以考虑先递归计算前一半，然后计算前一半对后一半的贡献，然后去递归后一半。

void CDQ(Poly& F, Poly& G, int l, int r) {
    if (l == r) {
        if (!l) F[l] = 1;
        return;
    }
    int mid = (l + r) >> 1;
    CDQ(F, G, l, mid);
    Poly A(F.begin() + l, F.begin() + mid + 1);
    Poly B(G.begin(), G.begin() + r - l + 1);
    A = A * B;
    for (int i = mid + 1; i <= r; ++i) 
        F[i] = add(F[i], A[i - l]);
    CDQ(F, G, mid + 1, r);
}

17. 完结撒花

大概是完了。

给个代码实现。

用 std::vector<int> 实现多项式全家桶其实并没有特别慢，跑的还是挺快的。而且代码看上去很舒服，不会特别乱糟糟的一大片。

小范围内暴力卷积会比 (FFT) 快不少。小范围大概就是 (2000) 左右这个样子。

#include <bits/stdc++.h>

inline int read() {
  int ret, cc, sign = 1;
  while (!isdigit(cc = getchar()))
    sign = cc == '-' ? -1 : sign;
  ret = cc - 48;
  while (isdigit(cc = getchar()))
    ret = cc - 48 + ret * 10;
  return ret;
}

const int MOD = 998244353;
const int G = 3;
const int MAXN = 600010;
typedef std::vector<int> Poly;

inline int add(int a, int b) { return (a += b) >= MOD ? a - MOD : a; }
inline int sub(int a, int b) { return (a -= b) <    0 ? a + MOD : a; }
inline int mul(int a, int b) { return 1ll * a * b % MOD; }
inline int qpow(int a, int p) {
  int ret = 1;
  for (p += (p < 0) * (MOD - 1); p; p >>= 1, a = mul(a, a))
    if (p & 1) ret = mul(ret, a);
  return ret;
}

inline Poly operator * (const Poly&, const Poly&);
inline Poly operator + (const Poly&, const Poly&);
inline Poly operator - (const Poly&, const Poly&);
inline Poly Inverse(const Poly&);
inline Poly Derivate(const Poly&);
inline Poly Integral(const Poly&);
inline Poly Ln(const Poly&);
inline Poly Exp(const Poly&);
inline Poly Pow(const Poly&, int);
inline Poly Sqrt(const Poly&);
inline Poly Multieval(const Poly&, const Poly&);
inline Poly Fastinsv(const Poly&, const Poly&);
inline std::pair<Poly, Poly> Div(const Poly&, const Poly&);
inline void Read(Poly&);
inline void Print(const Poly&);

int inv[MAXN];
int fac[MAXN];

int main() {

}

inline void Read(Poly& A) {
  std::generate(A.begin(), A.end(), read);
}

inline void Print(const Poly& A) {
  for (auto x : A) printf("%d ", x);
  puts("");
}

int rev[MAXN];
int W[MAXN];

inline int getrev(int n) {
  int len = 1, cnt = 0;
  while (len < n) len <<= 1, ++cnt;
  for (int i = 0; i < len; ++i)
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (cnt - 1));
  return len;
}

inline void NTT(Poly& a, int n, int opt) {
  a.resize(n);
  for (int i = 0; i < n; ++i)
    if (i < rev[i])
      std::swap(a[i], a[rev[i]]);
  for (int i = (W[0] = 1); i < n; i <<= 1) {
    int wn = qpow(G, opt * (MOD - 1) / (i << 1));
    for (int k = i - 2; k >= 0; k -= 2)
      W[k + 1] = mul(W[k] = W[k >> 1], wn);
    for (int j = 0, p = i << 1; j < n; j += p) {
      for (int k = 0; k < i; ++k) {
        int x = a[j + k], y = mul(W[k], a[j + k + i]);
        a[j + k] = add(x, y), a[j + k + i] = sub(x, y);
      }
    }
  }
  if (opt == -1)
    for (int i = 0, r = qpow(n, MOD - 2); i < n; ++i)
      a[i] = mul(a[i], r);
}

inline Poly operator * (const Poly& lhs, const Poly& rhs) {
  if (lhs.size() * rhs.size() <= 2000) {
    Poly A(lhs.size() + rhs.size() - 1);
    for (size_t i = 0; i < lhs.size(); ++i)
      for (size_t j = 0; j < rhs.size(); ++j)
        A[i + j] = add(A[i + j], mul(lhs[i], rhs[j]));
    return A;
  }
  Poly A(lhs), B(rhs);
  int len = A.size() + B.size() - 1, bln = getrev(len);
  NTT(A, bln, 1), NTT(B, bln, 1);
  for (int i = 0; i < bln; ++i)
    A[i] = mul(A[i], B[i]);
  NTT(A, bln, -1), A.resize(len);
  return A;
}

inline Poly operator + (const Poly& lhs, const Poly& rhs) {
  Poly C(std::max(lhs.size(), rhs.size()));
  for (size_t i = 0; i < C.size(); ++i)
    C[i] = add(i < lhs.size() ? lhs[i] : 0, i < rhs.size() ? rhs[i] : 0);
  return C;
}

inline Poly operator - (const Poly& lhs, const Poly& rhs) {
  Poly C(std::max(lhs.size(), rhs.size()));
  for (size_t i = 0; i < C.size(); ++i)
    C[i] = sub(i < lhs.size() ? lhs[i] : 0, i < rhs.size() ? rhs[i] : 0);
  return C;
}

inline Poly Inverse(const Poly& A) {
  Poly B(1, qpow(A[0], MOD - 2));
  int n = A.size() << 1;
  for (int i = 2; i < n; i <<= 1) {
    Poly C(A); C.resize(i);
    int len = getrev(i << 1);
    NTT(B, len, 1), NTT(C, len, 1);
    for (int j = 0; j < len; ++j)
      B[j] = mul(B[j], sub(2, mul(B[j], C[j])));
    NTT(B, len, -1), B.resize(i);
  }
  B.resize(A.size());
  return B;
}

inline std::vector<int> getinv() {
  std::vector<int> inv(MAXN);
  inv[1] = 1;
  for (int i = 2; i < MAXN; ++i)
    inv[i] = mul(MOD - MOD / i, inv[MOD % i]);
  return inv;
}

std::vector<int> Inv = getinv();

inline Poly Derivate(const Poly& A) {
  Poly C(A.size() - 1);
  for (size_t i = 0; i < C.size(); ++i)
    C[i] = mul(i + 1, A[i + 1]);
  return C;
}

inline Poly Integral(const Poly& A) {
  Poly C(A.size() + 1);
  for (size_t i = 1; i < C.size(); ++i)
    C[i] = mul(Inv[i], A[i - 1]);
  return C;
}

inline Poly Ln(const Poly& A) {
  Poly C = Integral(Derivate(A) * Inverse(A));
  C.resize(A.size());
  return C;
}

inline Poly Exp(const Poly& A) {
  Poly B(1, 1);
  int n = A.size() << 1;
  for (int i = 2; i < n; i <<= 1) {
    Poly C(A);
    C.resize(i), B.resize(i);
    B = B * (Poly(1, 1) - Ln(B) + C);
  }
  B.resize(A.size());
  return B;
}

inline Poly Pow(const Poly& A, int k) {
  Poly C(Ln(A));
  for (size_t i = 0; i < C.size(); ++i)
    C[i] = mul(C[i], k);
  return Exp(C);
}

inline Poly Sqrt(const Poly& A) {
  Poly C(A);
  int c = A[0], ic = qpow(c, MOD - 2);
  for (size_t i = 0; i < C.size(); ++i)
    C[i] = mul(C[i], ic);
  c = sqrt(c), C = Pow(C, Inv[2]);
  for (size_t i = 0; i < C.size(); ++i)
    C[i] = mul(C[i], c);
  return C;
}

inline std::pair<Poly, Poly> Div(const Poly& lhs, const Poly& rhs) {
  if (rhs.size() > lhs.size()) return {Poly(1, 0), lhs};
  Poly A(lhs), B(rhs);
  int len = A.size() - B.size() + 1;
  std::reverse(A.begin(), A.end());
  std::reverse(B.begin(), B.end());
  A.resize(len), B.resize(len);
  A = A * Inverse(B);
  A.resize(len);
  std::reverse(A.begin(), A.end());
  B = lhs - rhs * A;
  B.resize(rhs.size() - 1);
  return {A, B};
}

Poly mem[MAXN];
inline void Divide(int u, int l, int r, const Poly& A) {
  if (l == r) {
    mem[u] = Poly({MOD - A[l], 1});
    return;
  }
  int mid = (l + r) >> 1;
  Divide(u << 1, l, mid, A), Divide(u << 1 | 1, mid + 1, r, A);
  mem[u] = mem[u << 1] * mem[u << 1 | 1];
}

inline void Eval_recur(int u, int l, int r, const Poly& A, Poly& res) {
  if (l == r) return res.push_back(A[0]);
  int mid = (l + r) >> 1;
  Eval_recur(u << 1, l, mid, Div(A, mem[u << 1]).second, res);
  Eval_recur(u << 1 | 1, mid + 1, r, Div(A, mem[u << 1 | 1]).second, res);
}

inline Poly Multieval(const Poly& A, const Poly& B) {
  Divide(1, 0, B.size() - 1, B);
  Poly res;
  Eval_recur(1, 0, B.size() - 1, Div(A, mem[1]).second, res);
  return res;
}


inline Poly Insv_recur(int u, int l, int r, const Poly& A) {
  if (l == r) return Poly({A[l]});
  int mid = (l + r) >> 1;
  Poly a = Insv_recur(u << 1, l, mid, A) * mem[u << 1 | 1];
  Poly b = Insv_recur(u << 1 | 1, mid + 1, r, A) * mem[u << 1];
  return a + b;
}

inline Poly Fastinsv(const Poly& X, const Poly& Y) {
  Divide(1, 0, X.size() - 1, X);
  Poly A = Derivate(mem[1]), B;
  Eval_recur(1, 0, X.size() - 1, A, B);
  for (size_t i = 0; i < X.size(); ++i)
    B[i] = mul(qpow(B[i], MOD - 2), Y[i]);
  return Insv_recur(1, 0, X.size() - 1, B);
}