多项式全家桶

下面的所有知识均属多项式范畴。

---

前言

本文略去多项式的定义，毕竟初中数学都讲过。

本文会出现一些可能没有见过的表达：

• 对于一个关于 $x$ 的多项式 $f(x)$，有时也称其最高次项次数为 $f(x)$ 的度数。
• 多项式同样定义了四则运算，包括常数与多项式的四则运算和多项式之间的四则运算，同时也衍生出了同余符号 $\equiv$、对数函数 $\ln (f(x))$，指数函数 $\exp (f(x))$ 等等。
• 下文一般采用 $f(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$ 来表示多项式 $f(x)$ 的所有次数小于 $n$ 的项之和。
• 下文有时会采用 $[x^k]f(x)$ 表示多项式次数为 $k$ 的项。

注意：本文的知识环环相扣，请勿乱序阅读。

---

拉格朗日插值

插值，通俗地讲，就是指通过多项式在某些点的值求出整个多项式的表达式。

而拉格朗日插值就是计算插值的一种方法，它支持通过 $n+1$ 个点插出一个度数为 $n$ 的多项式。

按惯例，献上模板题。

首先对于多项式 $f(x)$，必然满足 $f(x)\equiv f(a)(\mathrm{m}\mathrm{o}\mathrm{d}(x-a))$。

证明：

转化为 $f(x)-f(a)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}(x-a))$。

设左侧的多项式为 $g(x)$，显然其常数项为 $0$。

又因为其第 $i$ 项为 $p_i(x^i-a^i)$（$p_i$ 为系数），并且 $x^i-a^i$ 必定可以分解为 $(x-a)\cdot q(x)$ 的形式（$q(x)$ 为另一个多项式），所以得证。

所以我们可以得到一个关于 $f(x)$ 线性同余方程组：
$$\{\begin{array}{c}f(x)\equiv y_1(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_1))\\ f(x)\equiv y_2(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_2))\\ f(x)\equiv y_3(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_3))\\ ⋮\\ f(x)\equiv y_n(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_n))\end{array}$$
这里可以使用中国剩余定理解方程组（这里就没有必要关注它们是否互质了）。

我们可以设 $M={\prod }_{i=1}^n(x-x_i)$，$M_i=\frac{M}{(x-x_i)}={\prod }_{j\ne i}(x-x_j)$。

则 $M_i$ 在模 $x-x_i$ 意义下的逆元为 ${\prod }_{j\ne i}\frac{1}{x_i-x_j}$。

（因为 $M_i{M}_i^{-1}={\prod }_{j\ne i}\frac{x-x_j}{x_i-x_j}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_i))$，则对于每一个 $j\ne i$ 都存在 $ x - x_j \equiv x_i - x_j~(\bmod (x - x_i))$，移项可得 $x-x_i\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_i))$，这显然是成立的。）

所以有 $f(x)={\sum }_{i=1}^n{y}_i{\prod }_{j\ne i}\frac{x-x_j}{x_i-x_j}(\mathrm{m}\mathrm{o}\mathrm{d}M)$，它在模意义下是唯一的。

（实际上有另外的证明方式证明了其在通常意义下也是唯一的，这里略去。）

本题中，我们只需要代入 $k$ 求解即可，时间复杂度为 $\mathcal{O}(n^2)$。不过更多时候，我们需要求出其所有系数，也可以做到 $\mathcal{O}(n^2)$，甚至还存在更低的复杂度，这里就不展开介绍了。

代码就不贴了。

---

多项式乘法

多项式乘法的定义，这里就不必多言了，初中数学都讲过。

代数形式：

设有多项式 $f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}_n$，$g(x)=b_0+b_{1}x+b_2{x}^2+\cdots +b_m{x}_m$。

则它们的乘法积为 $t(x)=c_0+c_{1}x+c_2{x}^2+\cdots +c_{n+m}{x}^{n+m}$。

其中，对于 $t(x)$ 的系数 $c_k$ 有公式：$c_k={\sum }_{i=0}^k{a}_i{b}_{k-i}$。

可以发现这个公式的形式十分类似卷积。但请注意：它们不是同一个东西，不要弄混淆了。

事实上，多项式乘法就是用于求卷积的。如果您对卷积感兴趣，建议百度。

---

现在，给出两个多项式 $f(x)$ 和 $g(x)$ 的系数，求出它们的乘法积 $t(x)$。$f(x)$ 和 $g(x)$ 的次数小于 $10^5$。

根据上面的公式，我们可以得到一个 $\mathcal{O}(n^2)$ 的优秀算法。

我们考虑换一种思想计算这个乘法积。

众所周知，有这样一条著名的定理：任意 $k$ 个点可以唯一确定一个至多为 $k-1$ 次的多项式。

所以我们得到了一个新的思路：先求出足够多的 $f(x)$ 和 $g(x)$ 的点值，将它们用 $\mathcal{O}(n)$ 的复杂度乘起来（横坐标不变，纵坐标相乘），再将这些点转化为一个多项式（上面提到了这样的过程叫做插值）。

这种方法叫做离散傅里叶变换（即著名的 $\mathrm{DFT}$）。

问题是，如果钦定一些 $x$ 去求对应的 $f(x)$ 的值，时间复杂度还是 $\mathcal{O}(n^2)$ 的。

所以，下面介绍快速傅里叶变换，同时会介绍快速数论变换（即著名的 $\mathrm{FFT}$ 和 $\mathrm{NTT}$）。

---

前置知识：单位根，原根。

还记得二次函数 $y=x^2$ 的性质吗？

它关于 $y$ 轴对称，所以我们如果代入 $x=k$，同样可以得到 $x=-k$ 时 $y$ 的值。

这启发我们，可以选取一些特殊的 $x$，使得我们可以在求出一些点值的同时，快速地求出另外的一些点值。

这个时候，单位根就派上用场了。

先讲一些使得单位根能够应用于此的性质：

1. 记 ${\omega }_n$ 为 $n$ 次单位根，即满足 ${\omega }_n^k=e^{\frac{2k\pi i}{n}}$，而 $e^{\frac{2k\pi i}{n}}=\cos \frac{2k\pi }{n}+i\sin \frac{2k\pi }{n}$，$i$ 为虚数单位。

（注意 $n$ 次单位根的 $n$ 一般取 $2$ 的正整数次幂，而下文皆如此。）

这个计算公式不仅可以让我们计算出单位根的具体数值，并且可以得出以下性质：

• ${\omega }_n^n=\cos 2\pi +i\sin 2\pi =1={\omega }_n^0$。

这条性质说明 $n$ 次单位根存在长度为 $n$ 的指数循环节；

• ${\omega }_n^{\frac{n}{2}}=\cos \pi +i\sin \pi =-1$。

这条性质间接地证明了 $n$ 次单位根的最小指数循环节长度为 $n$，同时为下文折半引理的证明做了铺垫。

2. 单位根消去引理：${\omega }_{2n}^{2k}={\omega }_n^k$。

根据公式有 ${\omega }_{2n}^{2k}=e^{\frac{4k\pi i}{2n}}=e^{\frac{2k\pi i}{n}}={\omega }_n^k$，得证。

这条性质可以做为沟通不同次单位根的桥梁，同时也证明了 $n$ 次单位根的值域是 $2n$ 次单位根值域的真子集。

3. 单位根折半引理：${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$。

根据上面的性质有 ${\omega }_n^{k+\frac{n}{2}}={\omega }_n^k\cdot {\omega }_n^{\frac{n}{2}}=-{\omega }_n^k$，得证。

这条性质可以做为沟通同次不同指数的单位根的桥梁，为 $\mathrm{DFT}$ 的分治优化创造了可能。

而我们该如何运用单位根的这些性质呢？

我们将 $f(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3$ 做一个小小的改动：
$$\begin{array}{rl}f(x)& =a_0+a_{1}x+a_2{x}^2+a_3{x}^3\\ & =a_0+a_2{x}^2+a_{1}x+a_3{x}^3\\ & =(a_0+a_2{x}^2)+x(a_1+a_3{x}^2)\\ & =f_1(x^2)+x{f}_2(x^2)\end{array}$$
可以发现，这样的变换，会使得原多项式被拆分为两个规模缩小了一半多项式。

所以，我们又有了一个新的算法模型：

我们考虑将 $f(x)$ 拆分为两个规模缩小了一半的多项式 $f_1(x),f_2(x)$，然后对它们进行分治，接着，我们可以用 $f_1(x),f_2(x)$ 在 $x=k^2$ 处的点值计算出 $f(x)$ 在 $x=k$ 处的点值。如果分治时每一层的复杂度都能够做到 $\mathcal{O}(n)$，那么总的时间复杂度就是 $\mathcal{O}(n{\log }_{2}n)$。

接下来，我们考虑利用单位根的性质，将每一层转化点值的操作复杂度变为 $\mathcal{O}(n)$。

如果我们将单位根 ${\omega }_n^k$ 和 ${\omega }_n^{k+\frac{n}{2}}$ 分别代入多项式 $f(x)$：
$$\begin{gathered} f({\omega }_n^k)=f_1({\omega }_n^{2k})+{\omega }_n^k{f}_2({\omega }_n^{2k})=f_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{f}_2({\omega }_{\frac{n}{2}}^k) \\ f({\omega }_n^{k+\frac{n}{2}})=f_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{f}_2({\omega }_n^{2k+n})=f_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{f}_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$
可以发现，只要我们求出了 $f_1(x),f_2(x)$ 在 $x={\omega }_{\frac{n}{2}}^k$ 处的点值，我们就可以直接得出 $f(x)$ 在 $x={\omega }_n^k,x={\omega }_n^{k+\frac{n}{2}}$ 两处的点值！

具体来说，我们只要在所有长度为 $n$ 的多项式中求出在所有 $n$ 次单位根处的点值即可。

这样，做到 $\mathcal{O}(n)$ 的复杂度已经易如反掌了。

我们再来理一遍这个算法吧：（同时加入了一些很多细节）

1. 我们找到一个 $p$，使得 $p$ 是 $2$ 的正整数次幂且 $p\ge n+m+1$。

这一步操作之前没有提到，但是很有必要，它为分治算法的进行提供了方便，这样我们就无需考虑多项式的度数；同时，这一步操作不会影响到点值的大小，也不会影响到最终插出的多项式的度数 （想一想，为什么）。

代码：（$q$ 的作用下文会提及）

int p = 1, q = -1;
while(p <= n + m) ++q, p <<= 1;

2. 接着，我们对两个多项式 $f(x)$ 和 $g(x)$ 分别跑 $\mathrm{FFT}$。在 $\mathrm{FFT}$ 运行过程中，如果当前多项式长度为 $1$，直接返回系数；否则，将当前的多项式 $f(x)$ 按照系数底数的奇偶性将项分类（这里请读者自己想一想为什么不可以直接对半分类），将其拆分为两个多项式，然后分治。

但是，如果这一步直接采用递归，时间似乎会有些吃不消。

我们不妨来探讨一下系数究竟是怎样分配的：（以 $n=8$ 为例）
{ a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 } ↓ { a 0 , a 2 , a 4 , a 6 }   { a 1 , a 3 , a 5 , a 7 } ↓ { a 0 , a 4 }   { a 2 , a 6 }   { a 1 , a 5 }   { a 3 , a 7 } ↓ { a 0 }   { a 4 }   { a 2 }   { a 6 }   { a 1 }   { a 5 }   { a 3 }   { a 7 } \{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\} \\ \downarrow \\ \{a_0,a_2,a_4,a_6\}~\{a_1,a_3,a_5,a_7\} \\ \downarrow \\ \{a_0, a_4\}~\{a_2, a_6\}~\{a_1, a_5\}~\{a_3, a_7\} \\ \downarrow \\ \{a_0\}~\{a_4\}~\{a_2\}~\{a_6\}~\{a_1\}~\{a_5\}~\{a_3\}~\{a_7\} {a0​,a1​,a2​,a3​,a4​,a5​,a6​,a7​}↓{a0​,a2​,a4​,a6​} {a1​,a3​,a5​,a7​}↓{a0​,a4​} {a2​,a6​} {a1​,a5​} {a3​,a7​}↓{a0​} {a4​} {a2​} {a6​} {a1​} {a5​} {a3​} {a7​}
把每个系数的底数分离：

$$0,4,2,6,1,5,3,7$$

转换成二进制：

$$000,100,010,110,001,101,011,111$$

我们再把它们翻转过来：
$$000,001,010,011,100,101,110,111$$

这不就是从小到大排序的吗？

所以我们就有了如下代码：（设 $r_i$ 为分配完系数后第 $i$ 个位置上的系数）

for(register int i = 0; i < p; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << q);

解释一下：第 $i$ 个位置的 $r$，不就等于第 $\lfloor \frac{i}{2}\rfloor$ 个位置的 $r$ 右移一位吗？

如果 $i$ 是奇数，我们再在最高位补上一个 $1$ 即可（$q$ 的作用就在这里）。

以上过程称为 $\mathrm{Rader}$ 排序。

3. 用拆分出的两个多项式的点值，求出当前多项式的点值。

代码：（$\mathrm{Complex}$ 为复数类，定义了构造函数，重定义了乘法）

inline void FFT(Complex *f, int n)
{
	for(register int i = 0; i < n; i++) if(i < rev[i]) swap(f[i], f[rev[i]]);
	for(register int len = 2; len <= n; len <<= 1)
	{
	   	int k = len >> 1;
		Complex _w(cos(PI / k), t * sin(PI / k));
		for(register int i = 0; i < n; i += len)
		{
			Complex w(1.0, 0.0);
			for(register int j = 0; j < k; j++)
			{
				Complex u = f[i + j], v = w * f[i + j + k];
				f[i + j] = u + v;
				f[i + j + k] = u - v;
           	     w = w * _w;
			}
		}
	}
	return;
}

4. 接着，我们将 $f(x)$ 和 $g(x)$ 对应的点值相乘，时间复杂度为 $\mathcal{O}(n)$。

你以为结束了？并没有，我们还有第五步：

5. 最后，我们将多项式插值出每一项的系数，算法结束。我们称这一步叫 $\mathrm{IFFT}$。

“算法结束” 四个字就在这里，近在咫尺，却又远在天边。

如果用拉格朗日插值计算，时间复杂度还是 $\mathcal{O}(n^2)$，功亏一篑。

难道前面的分析过程都只能打水漂？

不不不，我们冷静分析一波。

看看下面这张图：

多项式乘法的过程实际上类似上图中的矩阵乘法。

那么我们想得到后面的系数矩阵，不就只要把中间的这个矩阵求逆就行了吗！

而这个逆矩阵非常特殊：

它等于原矩阵的每个值求倒数，再乘上一个 $\frac{1}{n}$。

（然而笔者线性代数太菜了，无法解释其原因。）

既然无法解释，那就考虑推式子证明吧：

设 $y_k$ 是 $f(x),g(x)$ 在傅里叶变换后，在 $x={\omega }_p^k$ 处点值相乘的结果。

我们设多项式 $s(x)=y_0+y_{1}x+y_2{x}^2+\cdots +y_{p-1}{x}^{p-1}$ 在 $x={\omega }_n^0,{\omega }_n^{-1},\cdots \cdots ,{\omega }_n^{-(p-1)}$ 处的点值分别为 $u_{0\sim p-1}$，即满足：
$$u_k=\sum _{i=0}^{p-1}{y}_i({\omega }_p^{-k}{)}^i$$
现在我们要用它证明等式 $c_k=\frac{u_k}{p}$ 成立，其中 $c_k$ 指最终 $f(x)$ 和 $g(x)$ 的乘法积的系数。

我们直接莽式子：
$$\begin{array}{rl}{u}_k& =\sum _{i=0}^{p-1}{y}_i({\omega }_p^{-k}{)}^i\\ & =\sum _{i=0}^{p-1}(\sum _{j=0}^{p-1}{c}_j({\omega }_p^i{)}^j)({\omega }_p^{-k}{)}^i\\ & =\sum _{i=0}^{p-1}\sum _{j=0}^{p-1}{c}_j({\omega }_p^{j-k}{)}^i\\ & =\sum _{j=0}^{p-1}{c}_j\sum _{i=0}^{p-1}({\omega }_p^{j-k}{)}^i\end{array}$$
然后对 ${\sum }_{i=0}^{p-1}({\omega }_p^{j-k}{)}^i$ 这一部分分讨，设 $v=j-k$，则该式等于 ${\sum }_{i=0}^{p-1}({\omega }_p^v{)}^i$。

在这个式子的意义下，这就是一个等比数列求和。

如果公比 ${\omega }_p^v=1$，即 $v=0$，则此时原式等于 $p$；

否则，套用等比数列公式，原式等于：
$$\begin{array}{rl}& \frac{({\omega }^v{)}^p-1}{{\omega }_p^v-1}\\ =& \frac{({\omega }^p{)}^v-1}{{\omega }_p^v-1}\\ =& \frac{1-1}{{\omega }_p^v-1}\\ =& 0\end{array}$$
故当且仅当 $v=j-k=0$，即 $j=k$ 时，后面的这一部分才会产生贡献，则有 $u_k=p{c}_k$，即 $c_k=\frac{1}{p}{u}_k$。

得证。

然后可以发现，$\mathrm{FFT}$ 和 $\mathrm{IFFT}$ 的代码可以合为一份：

inline void FFT(Complex *f, int n, double t = 1.0)
{
	for(register int i = 0; i < n; i++) if(i < rev[i]) swap(f[i], f[rev[i]]);
	for(register int len = 2; len <= n; len <<= 1)
	{
		int k = len >> 1;
		Complex _w(cos(PI / k), t * sin(PI / k));
		for(register int i = 0; i < n; i += len)
		{
			Complex w(1.0, 0.0);
			for(register int j = 0; j < k; j++)
			{
				Complex u = f[i + j], v = w * f[i + j + k];
				f[i + j] = u + v;
				f[i + j + k] = u - v;
               	w = w * _w;
			}
		}
	}
   	if(t == -1) for(register int i = 0; i < n; i++) f[i].r /= n, f[i].i /= n;
	return;
}

当运行 $\mathrm{FFT}$ 时，我们令 $t=1$ 即可；当运行 $\mathrm{IFFT}$ 时，我们令 $t=-1$ 即可。

（注意：由于 $\cos$ 和 $\sin$ 的精度本身就十分地不稳定，加上复数乘法的运算次数较多，我们最后计算答案的时候需要处理这个误差。）

这样，算法真正大功告成啦！

---

优化：

可以发现，对于一个复数 $a+bi$，有 $(a+bi{)}^2=(a^2-b^2)+2abi$。

如果我们把 $g(x)$ 放到 $f(x)$ 的虚部上，对 $f(x)$ 跑一次 $\mathrm{FFT}$，然后将 $f(x)$ 的每个点值变为其平方。那么，虚部除以 $2$ 不就是我们想要的值吗？

这样，我们就可以只跑一遍 $\mathrm{FFT}$，一遍 $\mathrm{IFFT}$ 即可，常数优化至原来的 $\frac{2}{3}$。

---

然而，$\mathrm{FFT}$ 的精度仍然是一个令人头疼的问题。我们想考虑一种方式，使得整个运算建立在整数域上。

整数，想到数论，有一个东西恰好符合其所有要求——原根。

关于原根的定义和性质，可以参考我的博客整除和同余相关理论。

可以发现，使得单位根可用于 $\mathrm{FFT}$ 的性质，原根同样满足：

1. $n$ 次原根计算公式为 $g_n^k=g^{\frac{k(p-1)}{n}}\mathrm{m}\mathrm{o}\mathrm{d}p$，其中 $p$ 为模数，$g$ 特指模 $p$ 的最小原根。

这里的 $p$ 一般取 $998244353$，它的最小原根等于 $3$，并且有 $998244353=2^{23}\times 7\times 17+1$，其中 $2^{23}=8388608$，其大小完全可以满足我们的需求。

同样，我们也可以得出 $g_n^n=g_n^0\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 和 $g_n^{\frac{n}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 以及它们的推论。

2. 原根同样存在消去引理和折半引理。

所以，在 $\mathrm{FFT}$ 中对多项式的变换，在 $\mathrm{NTT}$ 中同样适用。同时，也需要 $\mathrm{Rader}$ 排序。

同理，对于 $\mathrm{INTT}$，原根同样可以在原来的基础上求倒数再乘 $\frac{1}{n}$，只需要求逆元即可。

代码：

#define MOD 998244353
#define Min_G 3
#define Inv_Min_G 332748118

inline int Get_MOD(int k) { if(k >= MOD) k -= MOD; return k; }

inline int Pow(int a, int b)
{
	int res = 1;
	while(b)
	{
		if(b & 1) res = 1LL * res * a % MOD;
		a = 1LL * a * a % MOD;
		b >>= 1;
    }	
	return res;
}

inline void NTT(int *f, int n, bool Inverse = false)
{
	for(register int i = 0; i < n; i++) if(i < rev[i]) swap(f[i], f[rev[i]]);
	for(register int len = 2; len <= n; len <<= 1)
	{
		int k = (len >> 1), g = Pow(Inverse ? Inv_Min_G : Min_G, (MOD - 1) / len);
		for(register int i = 0; i < n; i += len)
		{
			int w = 1;
			for(register int j = 0; j < k; j++)
			{
				int u = f[i + j], v = 1LL * w * f[i + j + k] % MOD;
				f[i + j] = Get_MOD(u + v);
				f[i + j + k] = Get_MOD(u - v + MOD); 
				w = 1LL * w * g % MOD;
			}
		}
	}
	if(Inverse) 
	{
		int Inv_N = Pow(n, MOD - 2);
		for(register int i = 0; i < n; i++) f[i] = 1LL * f[i] * Inv_N % MOD;
	}
	return;
}

---

后记：

实际上，$\mathrm{FWT}$ 的思想和它们几乎是一模一样的，不过 $\mathrm{FWT}$ 没有涉及到较多数学知识，会相对简单一些。

强推神佬的视频：快速傅里叶变换($\mathrm{FFT}$)——有史以来最巧妙的算法？

强推神佬的博客：$\mathrm{FFT}$ $\mathrm{NTT}$

---

多项式全家桶

到这里，正文才刚刚开始。。。

首先我们要知道一点求导和积分的前置知识。

对于要用到的结论，这里就直接摆出来了：
$$\begin{array}{rl}& f(x)=e^x\Rightarrow f^{\prime }(x)=e^x\\ & f(x)=\ln x\Rightarrow f^{\prime }(x)=\frac{1}{x}\\ & f(x)=a{x}^b\Rightarrow f^{\prime }(x)=ab{x}^{b-1}\\ & (f(x)\pm g(x){)}^{\prime }=f^{\prime }(x)\pm g^{\prime }(x)\\ & (f(g(x){)}^{\prime }=f^{\prime }(g(x))\ast g^{\prime }(x)\end{array}$$
其中 ${}^{\prime }$ 表示求导。

所以我们可以推出多项式求导公式：
$$f(x)=\sum _{i=0}^n{a}_i{x}^i\Rightarrow f^{\prime }(x)=\sum _{i=0}^{n-1}(i+1)a_{i+1}{x}^i$$
当然，如果知道了多项式的导数，我们也可以用积分求出原多项式：
$$f^{\prime }(x)=\sum _{i=0}^n{a}_i{x}^i\Rightarrow f(x)=\sum _{i=1}^{n+1}\frac{a_{i-1}{x}^i}{i}$$
注意：多项式积分是不能得到正确的常数项的。不过大部分题目在数据上做了保证，可以直接认为其为 $0$。

这两个过程的代码写起来就很简单了，时间复杂度都是 $\mathcal{O}(n)$ 的：

求导：（把度数为 $n-1$ 的多项式 $f(x)$ 求导，结果为 $g(x)$）

inline void Get_Derivative(int *f, int *g, int n)
{
	for(register int i = 1; i < n; i++)
        g[i - 1] = 1LL * i * f[i] % MOD;
	return g[n - 1] = 0, void();
}

积分：（把多项式 $f(x)$ 求积分，结果为度数为 $n-1$ 的多项式 $g(x)$）

inline void Get_Integral(int *f, int *g, int n)
{
	for(register int i = n - 1; i >= 1; i--)
        g[i] = 1LL * f[i - 1] * Pow(i, MOD - 2) % MOD;
	return g[0] = 0, void();
}

在文章后面的部分，有很多地方会把一些式子通过求导变成方便计算的形式，然后再用积分求出答案。

---

以上结论已经能够满足大部分 $\mathrm{OI}$ 题目的需求了。

如果对证明感兴趣，或者想知道更加严谨的结论，可以查找更多资料学习。（显然笔者也只会结论。。。）

---

多项式乘法逆元

给出一个度数为 $n-1$ 的多项式 $f(x)$，求出一个度数同为 $n-1$ 多项式 $g(x)$，满足 $f(x)\ast g(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

其中 $\ast$ 表示多项式乘法。系数对 $998244353$ 取模。

---

先给出算法：倍增。

假如我们已经求出了一个多项式 $t$，满足 $f(x)\ast t(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

又因为 $f(x)\ast g(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\Rightarrow f(x)\ast g(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

所以 $f(x)\ast (g(x)-t(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

$f(x)$ 由于要求逆元，所以不可能为 $0$，那么有 $g(x)-t(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

平方后移项有：$g^2(x)\equiv 2\ast g(x)\ast t(x)-t^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

然后两边再乘上一个 $f(x)$：$f(x)\ast g^2(x)\equiv 2\ast f(x)\ast g(x)\ast t(x)-f(x)\ast t^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

根据前面的结论，可以化简为：$g(x)\equiv 2\ast t(x)-f(x)\ast t^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

我们知道 $f(x)$ 和 $t(x)$，那么自然也可以求出 $g(x)$ 了。

对于常数项，直接求出其逆元即可。

时间复杂度 $\mathcal{O}(n{\log }_{2}n)$。

贴上代码：（注意数组之间的覆盖与清零）

inline void Init(int n)
{
	p = 1, q = -1;
	while(p < n) p <<= 1, ++q;
	for(register int i = 0; i < p; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << q);
	return;
}

inline void Get_Operation_Sequence(int n, int *seq)
{
	seq[0] = 0;
	while(n != 1)
	{
		seq[++seq[0]] = n;
		n = (n + 1) >> 1;
	}
	seq[++seq[0]] = 1;
	return;
}

inline void Get_Inverse_Element(int *f, int *g, int n)
{
	Get_Operation_Sequence(n, seq[0]);
	for(register int k = seq[0][0]; k >= 1; k--)
	{
		int m = seq[0][k];
		if(m == 1)
		{
			g[0] = Pow(f[0], MOD - 2);
			continue;
		}
		Init(m << 1);
		memcpy(t, f, sizeof(int) * m);
		fill(t + m, t + p, 0);
		NTT(t, p), NTT(g, p);
		for(register int i = 0; i < p; i++) g[i] = 1LL * g[i] * Get_MOD(2LL - 1LL * t[i] * g[i] % MOD + MOD) % MOD;
		NTT(g, p, true), fill(g + m, g + p, 0);
		continue;
	}
}

---

多项式除法

（话说多项式除法的地位和多项式乘法比起来为何如此卑微。。。）

给出一个度数为 $n-1$ 的多项式 $f(x)$ 和一个度数为 $m-1$ 多项式 $g(x)$，请求出多项式 $q(x)$ 和 $r(x)$。

满足：

• $f(x)=g(x)\ast q(x)+r(x)$；
• $q(x)$ 度数为 $n-m$，$r(x)$ 度数小于 $m-1$；

所有系数对 $998244353$ 取模。

---

没记错的话，初中数学老师讲过一种 “大除法”。

不过本题是在模意义下计算的，所以说这种方法不是很适用。

接下来，我们需要考虑用先前学过的知识

那么这操作就很关键了：系数翻转。

（为了方便，下文设 $f(x)$ 度数为 $n$，设 $g(x)$ 度数为 $m$，其余的度数变化同理。）

我们设对 $f(x)$ 的第 $0\sim n$ 个系数翻转后为 $f^{r(n)}(x)$。

（注意下文不会刻意地区分 $f^{r(n)}(x)$ 中 $r(n)$ 和余数 $r(x)$ 的区别。）

那么可以发现有：
$$\begin{gathered} f^{r(n)}(x)=\sum _{i=0}^n{a}_{n-i}{x}^i \\ \Downarrow \\ {f}^{r(n)}(x)=x^{n}f(\frac{1}{x}) \end{gathered}$$
又因为 $f(\frac{1}{x})=g(\frac{1}{x})\ast q(\frac{1}{x})+r(\frac{1}{x})$，所以：
$$\begin{array}{rl}{f}^{r(n)}(x)& =x^n\ast (g(x)\ast q(x)+r(x))\\ & =x^n\ast (g(x)\ast q(x))+x^{n}r(x)\\ & =x^m\ast g(x)\ast x^{n-m}\ast q(x)+x^{n}r(x)\\ & =g^{r(m)}(x)\ast q^{r(n-m)}(x)+r^{r(n)}(x)\end{array}$$
然后，既然 $r(x)$ 的度数小于 $m$，那么翻转过后，$r^{r(n)}(x)$ 的前 $n-m$ 个系数就必然为 $0$。

所以，如果我们在 $\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n-m+1}$ 意义下做上述运算，不仅可以去掉余数 $r(x)$ 的影响，同时我们还可以发现，这样完美地保留了 $q(x)$。

那么有：
$$\begin{gathered} f^{r(n)}(x)\equiv g^{r(m)}(x)\ast q^{r(n-m)}(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n-m+1}) \\ \Downarrow \\ {q}^{r(n-m)}(x)\equiv \frac{g^{r(m)}(x)}{f^{r(n)}(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n-m+1}) \\ \Downarrow \\ q(x)=(\frac{g^{r(m)}(x)}{f^{r(n)}(x)}{)}^{r(n-m)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n-m+1}) \end{gathered}$$
至于余数，我们直接把 $g(x)$ 和 $q(x)$ 乘起来，再用 $f(x)$ 减去即可。

代码：

inline void Multiply(int *f, int n, int *g, int m, int *v)
{
    Init(n + m - 1);
	memcpy(h, f, sizeof(int) * n), fill(h + n, h + p, 0);
	memcpy(s, g, sizeof(int) * m), fill(s + m, s + p, 0);
	NTT(h, p), NTT(s, p);
	for(register int i = 0; i < p; i++) v[i] = 1LL * h[i] * s[i] % MOD;
	return NTT(v, p, true), fill(v + n + m - 1, v + p, 0), void();
}

inline void Divide(int *f, int n, int *g, int m, int *v)
{
	memcpy(s, g, sizeof(int) * m);
	reverse(s, s + m);
	Get_Inverse_Element(s, r, n - m + 1);
	memcpy(t, f, sizeof(int) * n);
	reverse(t, t + n);
	Multiply(t, n, r, n - m + 1, v);
	return reverse(v, v + n - m + 1);
}

inline void Get_Remainder(int *f, int n, int *g, int m, int *v, int *x)
{
	Multiply(g, m, v, n - m + 1, x);
	for(register int i = 0; i < m; i++) x[i] = Get_MOD(f[i] - x[i] + MOD);
	return;
}

---

多项式开根

给出一个度数为 $n-1$ 的多项式 $f(x)$，求出一个度数同为 $n-1$ 多项式 $g(x)$，满足 $g^2(x)\equiv f(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

保证有 $[x^0]f(x)=1$。系数对 $998244353$ 取模。

---

我们还是用倍增法来做。

临界条件：$n=1$ 时，$[x^0]g(x)=1$。

设已经求出了一个多项式 $t(x)$，满足 $t^2(x)\equiv f(x(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

同时设 $g(x{)}^2\equiv f(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

然后又是一波推式子：
$$\begin{array}{rl}& t^2(x)\equiv f(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })\\ & t^2(x)-f(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })\\ & (t^2(x)-f(x){)}^2\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ & (t^2(x)+f(x){)}^2\equiv 4\ast t^2(x)\ast f(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ & (\frac{t^2(x)+f(x)}{2\ast t(x)}{)}^2\equiv f(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ & \frac{t^2(x)+f(x)}{2\ast t(x)}\equiv g(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\\ & \frac{t(x)}{2}+\frac{f(x)}{2\ast t(x)}\equiv g(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)\end{array}$$
推到这里式子就好做了。时间复杂度 $\mathcal{O}(n{\log }_{2}n)$。

当然还有一种常数较小的写法：在维护 $t(x)$ 的同时维护 $\frac{1}{t(x)}$，因为开根的递归结构和求逆是一样的。

贴上代码：（没有写这种常数较小的写法）

inline void Get_Sqrt(int *f, int *g, int n)
{
	Get_Operation_Sequence(n, seq[1]);
	for(register int k = seq[1][0]; k >= 1; k--)
	{
		int m = seq[1][k];
		if(m == 1)
		{
			g[0] = Quadratic_Residue :: Solve(f[0]);
			continue;
		}
		Init(m << 1);
		memset(h, 0, sizeof(int) * p);
		Get_Inverse_Element(g, h, m);
		memcpy(s, f, sizeof(int) * m);
		fill(s + m, s + p, 0);
		NTT(g, p), NTT(h, p), NTT(s, p);
		int Inv2 = Pow(2, MOD - 2);
		for(register int i = 0; i < p; i++) g[i] = 1LL * Inv2 * Get_MOD(g[i] + 1LL * h[i] * s[i] % MOD) % MOD;
		NTT(g, p, true), fill(g + m, g + p, 0);
		continue;
	}
	return;
}

---

当然，还可以去做一下加强版，仅需要在求常数项时套上一个二次剩余即可。

写好二次剩余后，就只需要改 Get_Sqrt函数的第一行即可。

代码就不放出来了。

---

多项式对数函数

给出一个度数为 $n-1$ 的多项式 $f(x)$，求出一个度数同为 $n-1$ 多项式 $g(x)$，满足 $g(x)\equiv \ln (f(x))(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

保证有 $[x^0]f(x)=1$。系数对 $998244353$ 取模。

---

看到这里，大多数读者应该都是一脸懵：怎么对一个多项式求自然对数啊？

不要急，还记得前置知识吗？求导和积分。

前面还提到了一句话：把一些式子通过求导变成方便计算的形式，然后再用积分求出答案。

那么，做法已经呼之欲出了：

对 $g(x)$ 求导有：$g^{\prime }(x)\equiv {\ln }^{\prime }(f(x))\ast f^{\prime }(x)\equiv \frac{f^{\prime }(x)}{f(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

这不就是求一个导数和一个逆元吗？

最后积分回去就好啦。

对于常数项，显然有 $e^0=1$，即 $\ln 1=0$，所以保证了 $[x^0]f(x)=1$，否则是很难做的。

时间复杂度 $\mathcal{O}(n{\log }_{2}n)$。

贴上代码：

inline void Get_Ln(int *f, int *g, int n)
{
	Get_Inverse_Element(f, g, n);
	memset(t, 0, sizeof(int) * p);
	Get_Derivative(f, t, n);
	NTT(g, p), NTT(t, p);
	for(register int i = 0; i < p; i++) g[i] = 1LL * g[i] * t[i] % MOD;
	return NTT(g, p, true), Get_Integral(g, g, n), void();
}

---

多项式牛顿迭代

给出一个度数为 $n-1$ 的多项式 $f(x)$，求出一个度数同为 $n-1$ 多项式 $g(x)$，满足 $f(g(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

保证有 $[x^0]f(x)=0$。系数对 $998244353$ 取模。

---

首先需要知道前置知识：泰勒展开。

对于多项式 $f(x)$，如果我们已经知道了 $f(x_0)$ 的值，那么 $f(x)$ 的值可以用如下方法表示：
$$f(x)=\sum _{i=0}^{\infty }(\frac{f^{(i)}(x_0)}{i!}\ast (x-x_0{)}^i)=f(x_0)+\sum _{i=1}^{\infty }(\frac{f^{(i)}(x_0)}{i!}\ast (x-x_0{)}^i)$$
其中 $f^{(i)}(x_0)$ 指的是多项式 $f(x)$ 在 $x_0$ 处的 $i$ 阶导数，简单来讲就是对它求 $i$ 次导数得出的值。

这个公式的证明和其它用法大可不必知道，我们只会用它做一些最基本的操作，不会超出我们的认知范围。

---

还是考虑倍增。

对于 $n=1$ 的情况，此时我们一般会默认 $[x^0]g(x)=1$。

设求出了一个多项式 $t(x)$，满足 $f(t(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

那么我们要求多项式 $g(x)$，满足 $f(g(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

又可以得出：$f(g(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

如果您阅读完了上面的《拉格朗日插值》一节，就一定可以知道这个定理：

对于多项式 $f(x)$，必然满足 $f(x)\equiv f(a)(\mathrm{m}\mathrm{o}\mathrm{d}(x-a))$。

那么我们很快可以得出结论：$g(x)-t(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$。

而这个结论有什么作用呢？往下看。

我们要求的是 $g(x)$，满足 $f(g(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

那么我们在 $t(x)$ 处进行泰勒展开，得到：

$$f(g(x))\equiv f(t(x))+\sum _{i=1}^{\infty }(\frac{f^{(i)}(t(x))}{i!}\ast (g(x)-t(x){)}^i)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$

又由于上面的结论，我们可以发现，在 $(g(x)-t(x){)}^i$ 这个式子中，当 $i\ge 2$ 时，前 $n$ 位必然为 $0$。

所以整个式子就等于：
$$f(g(x))\equiv f(t(x))+f^{\prime }(t(x))\ast (g(x)-t(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
稍加推导一下就有：
$$g(x)=t(x)-\frac{f(t(x))}{f^{\prime }(t(x))}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
这样就可以倍增求出 $g(x)$ 了。时间复杂度 $\mathcal{O}(n{\log }_{2}n)$。

---

多项式指数函数

给出一个度数为 $n-1$ 的多项式 $f(x)$，求出一个度数同为 $n-1$ 多项式 $g(x)$，满足 $g(x)\equiv e^{f(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

保证有 $[x^0]f(x)=0$。系数对 $998244353$ 取模。

---

显然把 $e^{f(x)}$ 求导不会有任何变化。

所以，我们考虑用刚刚学到的牛顿迭代法做这道题。

我们构造出一个多项式（不管它的实际系数）$h(g(x))\equiv \ln (g(x))-f(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

那么把 $g(x)\equiv e^{f(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$ 代入就有：
$$h(g(x))\equiv \ln (g(x))-f(x)\equiv \ln (e^{f(x)})-f(x)\equiv f(x)-f(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
然后我们就可以套用上面的多项式牛顿迭代了：
$$g(x)\equiv t(x)-\frac{\ln (t(x))-f(x)}{\frac{1}{t(x)}}\equiv t(x)-t(x)\ast (\ln (t(x))-f(x))(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$$
注意一点：$f(x)$ 在求导的时候是被看作常数的。

对于常数项，可以参考上文《多项式对数函数》的证明。

时间复杂度 $\mathcal{O}(n{\log }_{2}n)$。

代码：

inline void Get_Exp(int *f, int *g, int n)
{
	Get_Operation_Sequence(n, seq[2]);
	for(register int k = seq[2][0]; k >= 1; k--)
	{
		int m = seq[2][k];
		if(m == 1)
		{
			g[0] = 1;
			continue;
		}
		Get_Ln(g, h, m);
		Init(m << 1);
		fill(h + m, h + p, 0);
		memcpy(s, f, sizeof(int) * m);
		fill(s + m, s + p, 0);
		NTT(g, p), NTT(h, p), NTT(s, p);
		for(register int i = 0; i < p; i++) g[i] = 1LL * g[i] * Get_MOD(Get_MOD(1LL - h[i] + MOD) + s[i]) % MOD;
		NTT(g, p, true), fill(g + m, g + p, 0), memset(h, 0, sizeof(int) * p), void();
	}
}

---

学完了这里，您可以尝试用其去推一推多项式求逆和多项式开根的式子，可以发现和原来推出来的式子是一样的。

---

多项式快速幂

给出一个度数为 $n-1$ 的多项式 $f(x)$ 和一个正整数 $k$，其中 $k\le 10^{10^5}$，求出一个度数同为 $n-1$ 多项式 $g(x)$，满足 $g(x)\equiv e^{f(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

保证有 $[x^0]f(x)=1$。系数对 $998244353$ 取模。

---

如果刷高中数学题刷得有感觉，应该可以一眼看出这里可以求 $\ln$ 后求 $\exp$。

所以，我们求出 $t(x)=\ln (f(x))$，那么有 $\ln (g(x))=k\ast t(x)$。

于是我们把 $t(x)$ 的每一项的系数乘上一个 $k$，最后对 $t(x)$ 求 $\exp$，得到的即为 $g(x)$。

时间复杂度 $\mathcal{O}(n{\log }_{2}n)$。

---

不过在贴代码前，有一个地方需要思考：$k$ 在读入时应该对什么取模？

可能会有人脱口而出是 $\phi (998244353)$。但是，需要郑重声明的是：欧拉定理在多项式的运算意义下是不成立的，正确的模数应当为 $998244353$ 本身（阅读完下文后，可以发现这也是规定了 $[x^0]f(x)=1$ 的原因）。

证明：

我们考虑证明一个更一般的情况：$f^p(x)\equiv f(x^p)(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

由二项式定理和 $Lucas$ 定理可以得出：$(a+b{)}^p\equiv a^p+b^p(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

那么由数学归纳法，我们设 $f(x)$ 度数为 $k$ 且 $f(x)=g(x)+a_k{x}^k$，同时上述结论对于度数为 $k-1$ 的多项式成立。

则有 $f^p(x)\equiv g^p(x)+a_k^p{x}^{kp}\equiv g(x^p)+a_k{x}^{kp}\equiv f(x^p)(\mathrm{m}\mathrm{o}\mathrm{d}p)$。（这里可以使用欧拉定理）

在这道题中，由于 $n<p$，所以在模 $p$ 意义下有 $f^p(x)\equiv f(x^p)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

由因为 $[x^0]f(x)=1$，所以在模 $p$ 意义下 $f^p(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

---

代码：

inline int Read_Get_MOD()
{
	int res = 0; char c = getchar();
	while(!isdigit(c)) c = getchar();
	while(isdigit(c))
	{
		res = (10LL * res + (c - '0')) % MOD; 
		c = getchar();
	}
	return res;
}

inline void Get_Quick_Pow(int *f, int *g, int n, int k)
{
	Get_Ln(f, r, n);
	for(register int i = 0; i < p; i++) r[i] = 1LL * r[i] * k % MOD;
	return Get_Exp(r, g, n);
}

---

加强版不保证 $[x^0]f(x)=1$？

我们如果把所有数除去常数项，就可以仿照上述方法了。当然，不要忘记在最后把这个除掉的数乘回去（注意这个数是整数，根据费马小定理，它的指数需要对 $\phi (998244353)$ 取模）。

不过，当然会存在常数项为 $0$ 的情况。这也很好解决，我们只要找到第一个 $a_t$（系数）不为 $0$ 的 $t$，就可以套用如下公式：
$$f^k(x)=x^{tk}\ast (\frac{f(x)}{x^t}{)}^k$$
从本质上讲，这个操作就是降次之后升次。

注意这个式子中的 $k$ 还是对 $998244353$ 取模。

贴上代码：（由于常数过大，不吸氧只有 $80$ 分，不过笔者实在是卡不动了。。。）

inline void Read_Get_MOD(int &k1, int &k2, int &k3)
{
	k1 = k2 = k3 = 0;
	char c = getchar();
	while(!isdigit(c)) c = getchar();
	while(isdigit(c))
	{
		k1 = (10LL * k1 + (c - '0')) % MOD; 
		k2 = (10LL * k2 + (c - '0')) % (MOD - 1);
		if(k3 < n) k3 = 10 * k3 + (c - '0'); 
		c = getchar();
	}
	return;
}

inline void Get_Quick_Pow(int *f, int *g, int n, int k1, int k2)
{
	int st = 0;
	while(st < n && f[st] == 0) ++st;
	if(st == n || 1LL * st * k1 > n)
	{
		for(register int i = 0; i < n; i++) putchar('0'), putchar(' ');
		putchar('\n');
		return;
	}
	int inv = Pow(f[st], MOD - 2), mul = Pow(f[st], k2);
	for(register int i = 0; i < n; i++) f[i] = 1LL * f[i + st] * inv % MOD;
	Get_Ln(f, r, n);
	for(register int i = 0; i < p; i++) r[i] = 1LL * r[i] * k1 % MOD;
	Get_Exp(r, g, n);
	int tk = st * k1;
	for(register int i = 0; i < tk; i++) putchar('0'), putchar(' ');
	for(register int i = tk; i < n; i++) printf("%lld ", 1LL * g[i - tk] * mul % MOD);
	putchar('\n');
	return;
}

inline void Run()
{
    int k1, k2, k3;
	scanf("%d", &n);
	Read_Get_MOD(k1, k2, k3);
	for(register int i = 0; i < n; i++) scanf("%d ", f + i);
	if(f[0] == 0 && k3 >= n)
	{
		for(register int i = 0; i < n; i++) putchar('0'), putchar(' ');
		putchar('\n');
		return;
	}
	Get_Quick_Pow(f, g, n, k1, k2);
	return;	
}

---

最后，附上这篇文章的总代码。

除快速幂的加强版需要吸氧之外，上文提到的所有模板题均可通过，仅需改写 Run()函数即可。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<cctype>
#include<ctime>
#include<algorithm>

using namespace std;

struct Complex
{
	double r, i;
	
	Complex(double _r = 0.0, double _i = 0.0) : r(_r), i(_i) {}
	
	Complex operator + (const Complex &t) const { return Complex(r + t.r, i + t.i); }
	Complex operator - (const Complex &t) const { return Complex(r - t.r, i - t.i); }
	Complex operator * (const Complex &t) const { return Complex(r * t.r - i * t.i, r * t.i + i * t.r); }
};

namespace Quadratic_Residue
{
	#define MOD 998244353
	
	int si;

	struct Complex_MOD
	{
		int a, b;

		Complex_MOD operator * (const Complex_MOD &t) const
		{
			Complex_MOD res;
			res.a = (1LL * a * t.a % MOD + 1LL * b * t.b % MOD * si % MOD) % MOD;
			res.b = (1LL * a * t.b % MOD + 1LL * b * t.a % MOD) % MOD;
			return res;
		}
	};

	inline int Pow_Real(int a, int b)
	{
		int res = 1;
		while(b)
		{
			if(b & 1) res = 1LL * res * a % MOD;
			a = 1LL * a * a % MOD;
			b >>= 1;
		}
		return res;
	}

	inline int Pow_Complex(Complex_MOD val, int t)
	{
		Complex_MOD res;
		res.a = 1, res.b = 0;
		while(t)
		{
			if(t & 1) res = res * val;
			val = val * val;
			t >>= 1;
		}
		return res.a;
	}

	inline int Solve(int n)
	{
		if(n == 0) return 0;
		srand((unsigned)(time(NULL)));
		int p1, p2;
		while(true)
		{
			p1 = 1LL * rand() * rand() % MOD;
			p2 = (1LL * p1 * p1 % MOD - n + MOD) % MOD;
			if(Pow_Real(p2, (MOD - 1) / 2) != 1) break;
		}
		si = p2;
		Complex_MOD val;
		val.a = p1, val.b = 1;
		int ans = Pow_Complex(val, (MOD + 1) / 2);
		if(MOD - ans < ans) ans = MOD - ans;
		return ans;
	}

	#undef MOD
}

namespace Polynomial
{
	#define MAXN 5000005
	#define MOD 998244353
	#define Min_G 3
	#define Inv_Min_G 332748118

	int n, f[MAXN], g[MAXN], t[MAXN], h[MAXN], s[MAXN], r[MAXN], u[MAXN];
	int p, q, rev[MAXN], seq[3][MAXN];
	const double PI = acos(-1.0);

	inline int Get_MOD(int k) { if(k >= MOD) k -= MOD; return k; }

	inline void Read_Get_MOD(int &k1, int &k2, int &k3)
	{
		k1 = k2 = k3 = 0;
		char c = getchar();
		while(!isdigit(c)) c = getchar();
		while(isdigit(c))
		{
			k1 = (10LL * k1 + (c - '0')) % MOD; 
			k2 = (10LL * k2 + (c - '0')) % (MOD - 1);
			if(k3 < n) k3 = 10 * k3 + (c - '0'); 
			c = getchar();
		}
		return;
	}

	inline int Pow(int a, int b)
	{
		int res = 1;
		while(b)
		{
			if(b & 1) res = 1LL * res * a % MOD;
			a = 1LL * a * a % MOD;
			b >>= 1;
		}
		return res;
	}

	inline void Init(int n)
	{
		p = 1, q = -1;
		while(p < n) p <<= 1, ++q;  
		for(register int i = 0; i < p; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << q);
		return;
	}

	inline void Get_Operation_Sequence(int n, int *seq)
	{
		seq[0] = 0;
		while(n != 1)
		{
			seq[++seq[0]] = n;
			n = (n + 1) >> 1;
		}
		seq[++seq[0]] = 1;
		return;
	}

	inline void Get_Derivative(int *f, int *g, int n)
	{
		for(register int i = 1; i < n; i++)
			g[i - 1] = 1LL * i * f[i] % MOD;
		return g[n - 1] = 0, void();
	}

	inline void Get_Integral(int *f, int *g, int n)
	{
		for(register int i = n - 1; i >= 1; i--)
			g[i] = 1LL * f[i - 1] * Pow(i, MOD - 2) % MOD;
		return g[0] = 0, void();
	}

	inline void FFT(Complex *f, int n, double t = 1.0)
	{
		for(register int i = 0; i < n; i++) if(i < rev[i]) swap(f[i], f[rev[i]]);
		for(register int len = 2; len <= n; len <<= 1)
		{
	   	    int k = len >> 1;
			Complex _w(cos(PI / k), t * sin(PI / k));
			for(register int i = 0; i < n; i += len)
			{
				Complex w(1.0, 0.0);
				for(register int j = 0; j < k; j++)
				{
					Complex u = f[i + j], v = w * f[i + j + k];
					f[i + j] = u + v;
					f[i + j + k] = u - v;
           	     	w = w * _w;
				}
			}
		}
  	  	if(t == -1) for(register int i = 0; i < n; i++) f[i].r /= n, f[i].i /= n;
		return;
	}

	inline void NTT(int *f, int n, bool Inverse = false)
	{
		for(register int i = 0; i < n; i++) if(i < rev[i]) swap(f[i], f[rev[i]]);
		for(register int len = 2; len <= n; len <<= 1)
		{
			int k = (len >> 1), g = Pow(Inverse ? Inv_Min_G : Min_G, (MOD - 1) / len);
			for(register int i = 0; i < n; i += len)
			{
				int w = 1;
				for(register int j = 0; j < k; j++)
				{
					int u = f[i + j], v = 1LL * w * f[i + j + k] % MOD;
					f[i + j] = Get_MOD(u + v);
					f[i + j + k] = Get_MOD(u - v + MOD); 
					w = 1LL * w * g % MOD;
				}
			}
		}
		if(Inverse) 
		{
			int Inv_N = Pow(n, MOD - 2);
			for(register int i = 0; i < n; i++) f[i] = 1LL * f[i] * Inv_N % MOD;
		}
		return;
	}

	inline void Get_Inverse_Element(int *f, int *g, int n)
	{
		Get_Operation_Sequence(n, seq[0]);
		for(register int k = seq[0][0]; k >= 1; k--)
		{
			int m = seq[0][k];
			if(m == 1)
			{
				g[0] = Pow(f[0], MOD - 2);
				continue;
			}
			Init(m << 1);
			memcpy(t, f, sizeof(int) * m);
			fill(t + m, t + p, 0);
			NTT(t, p), NTT(g, p);
			for(register int i = 0; i < p; i++) g[i] = 1LL * g[i] * Get_MOD(2LL - 1LL * t[i] * g[i] % MOD + MOD) % MOD;
			NTT(g, p, true), fill(g + m, g + p, 0);
			continue;
		}
	}
	
	inline void Multiply(int *f, int n, int *g, int m, int *v)
	{
		Init(n + m - 1);
		memcpy(h, f, sizeof(int) * n), fill(h + n, h + p, 0);
		memcpy(s, g, sizeof(int) * m), fill(s + m, s + p, 0);
		NTT(h, p), NTT(s, p);
		for(register int i = 0; i < p; i++) v[i] = 1LL * h[i] * s[i] % MOD;
		return NTT(v, p, true), fill(v + n + m - 1, v + p, 0), void();
	}

	inline void Divide(int *f, int n, int *g, int m, int *v)
	{
		memcpy(s, g, sizeof(int) * m);
		reverse(s, s + m);
		Get_Inverse_Element(s, r, n - m + 1);
		memcpy(t, f, sizeof(int) * n);
		reverse(t, t + n);
		Multiply(t, n, r, n - m + 1, v);
		return reverse(v, v + n - m + 1);
	}

	inline void Get_Remainder(int *f, int n, int *g, int m, int *v, int *x)
	{
		Multiply(g, m, v, n - m + 1, x);
		for(register int i = 0; i < m; i++) x[i] = Get_MOD(f[i] - x[i] + MOD);
		return;
	}

	inline void Get_Sqrt(int *f, int *g, int n)
	{
		Get_Operation_Sequence(n, seq[1]);
		for(register int k = seq[1][0]; k >= 1; k--)
		{
			int m = seq[1][k];
			if(m == 1)
			{
				g[0] = Quadratic_Residue :: Solve(f[0]);
				continue;
			}
			Init(m << 1);
			memset(h, 0, sizeof(int) * p);
			Get_Inverse_Element(g, h, m);
			memcpy(s, f, sizeof(int) * m);
			fill(s + m, s + p, 0);
			NTT(g, p), NTT(h, p), NTT(s, p);
			int Inv2 = Pow(2, MOD - 2);
			for(register int i = 0; i < p; i++) g[i] = 1LL * Inv2 * Get_MOD(g[i] + 1LL * h[i] * s[i] % MOD) % MOD;
			NTT(g, p, true), fill(g + m, g + p, 0);
			continue;
		}
		return;
	}

	inline void Get_Ln(int *f, int *g, int n)
	{
		Get_Inverse_Element(f, g, n);
		memset(t, 0, sizeof(int) * p);
		Get_Derivative(f, t, n);
		NTT(g, p), NTT(t, p);
		for(register int i = 0; i < p; i++) g[i] = 1LL * g[i] * t[i] % MOD;
		return NTT(g, p, true), Get_Integral(g, g, n), void();
	}

	inline void Get_Exp(int *f, int *g, int n)
	{
		Get_Operation_Sequence(n, seq[2]);
		for(register int k = seq[2][0]; k >= 1; k--)
		{
			int m = seq[2][k];
			if(m == 1)
			{
				g[0] = 1;
				continue;
			}
			Get_Ln(g, h, m);
			Init(m << 1);
			fill(h + m, h + p, 0);
			memcpy(s, f, sizeof(int) * m);
			fill(s + m, s + p, 0);
			NTT(g, p), NTT(h, p), NTT(s, p);
			for(register int i = 0; i < p; i++) g[i] = 1LL * g[i] * Get_MOD(Get_MOD(1LL - h[i] + MOD) + s[i]) % MOD;
			NTT(g, p, true), fill(g + m, g + p, 0), memset(h, 0, sizeof(int) * p), void();
		}
	}

	inline void Get_Quick_Pow(int *f, int *g, int n, int k1, int k2)
	{
		int st = 0;
		while(st < n && f[st] == 0) ++st;
		if(st == n || 1LL * st * k1 > n)
		{
			for(register int i = 0; i < n; i++) putchar('0'), putchar(' ');
			putchar('\n');
			return;
		}
		int inv = Pow(f[st], MOD - 2), mul = Pow(f[st], k2);
		for(register int i = 0; i < n; i++) f[i] = 1LL * f[i + st] * inv % MOD;
		Get_Ln(f, r, n);
		for(register int i = 0; i < p; i++) r[i] = 1LL * r[i] * k1 % MOD;
		Get_Exp(r, g, n);
		int tk = st * k1;
		for(register int i = 0; i < tk; i++) putchar('0'), putchar(' ');
		for(register int i = tk; i < n; i++) printf("%lld ", 1LL * g[i - tk] * mul % MOD);
		putchar('\n');
		return;
	}

	inline void Run()
	{
		/*
		 * You wanna do something here?
		 */
	}

	#undef MAXN
	#undef MOD
	#undef Min_G
	#undef Inv_Min_G
}

int main()
{
	Polynomial :: Run();
	return 0;
}

---

$ToBeContinued\cdots \cdots$