多项式操作

约定不等于承诺〃 2021-10-29 07:58 462阅读 0赞

*  多项式运算从入门到入坟
    
     *  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\$ 的时候，这个式子的值永远不会小于 $\\frac 1 2 at^2 $ 。而这一道永远迈不过去的，真实存在的 “坎”，就是极限。这就是微积分中最基础的思想——极限。

以上过程实际上就是积分的过程，用更正规的表达方式来表示，就是：  
\\\[ 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 泰勒展开和泰勒级数 ###

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

已知 \$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);
    }

转载于:https://www.cnblogs.com/chiaki-cage/p/11184972.html

[Link 1]: https://loj.ac/problem/2289