多项式运算JAVA 多项式运算是几年级

转载

mob64ca141a683a 2024-01-19 16:11:57

文章标签 多项式运算JAVA 多项式不定积分插值 文章分类 Java 后端开发

多项式全家桶

多项式运算从入门到入坟

多项式运算从入门到入坟

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 = ak\Delta t$

\[\begin {align} x=&\sum_{k = 1}^{N} v_k\Delta t \\ =& \sum_{k = 1}^{N} ak\Delta t^2 \\ =& a\Delta t^2\sum_{k = 1}^Nk \\ =& a\Delta t^2 \frac {N(N+1)}{2} \\ =& \frac 1 2 a\Delta t^2 \left(\left(\frac {t}{\Delta t}{}\right)^2 + \left(\frac {t}{\Delta t}{}\right)\right) \\ =& \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$

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

\[x = \int_0^tv\mathrm 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$

\[\begin {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'$

\[\begin {align} c' &= 0 \\ x^n {'} &= nx^{n-1} \\ a^x {'} &= a^x\ln 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$

\[\begin {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)}\right)' &= \frac {f'(x)g(x) - f(x)g'(x)} {g(x)^2} \\ (f(g(x)))' &= f'(g(x))g'(x) \end {align} \]

2. 6 常用积分公式

只需要知道两个。

\[\begin {align} \int x^n\mathrm dx &= \frac {1} {n+1}x^{n+1} + C \ (n \ne -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))$

\[\begin {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$

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

\[\begin {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^{k\ln A(x)} $ 。

10. 多项式开方

你可以选择牛顿迭代推式子，你也可以选择用上面的快速幂，也就是 $\text{Inv}2$

那么 $\text{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. 多项式带余除法

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

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

这里设

\[\begin{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)$

那我们尝试推一下式子。

\[\begin {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$

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

首先根据泰勒展开我们有

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

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

所以

\[\begin{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} \\ \]

然后

\[\begin{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(n\log^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}\neq{i}}{({x}-{x_j})}}{\prod_{{j}\neq{i}}{({x_i}-{x_j})}}{y_i}} \]

分母是个常数，单独和 $y$

\[F(x)=\sum_{i=1}^n{\frac{y_i}{\prod_{{j}\neq{i}}{({x_i}-{x_j})}}\prod_{{j}\neq{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}\neq{i}}{({x_i}-{x_j})}}$

后面的那个 $\prod_{{j}\neq{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)}\right)^n \]

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

\[[x^n]h(g(x))=\frac 1 n[x^{n-1}]h'(x)\left(\frac x {f(x)}\right)^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);
}

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。

上一篇：android蓝牙通信开发蓝牙通讯开发

下一篇：sql server 2012 sa无法登录 sql server2012用户sa登录失败

提问和评论都可以，用心的回复会被更多人看到评论

发布评论

相关文章

官方博客	全部文章	热门标签	班级博客
了解我们	网站地图	意见反馈

鸿蒙开发者社区	51CTO学堂
51CTO	软考资讯