多项式：从入门到全家桶

第一次用博客园写学习笔记，写的不好请见谅~

欢迎各位OIer在评论区说一下自己对这篇博客的建议！

有什么问题也欢迎在下方评论区中提出！顺便点赞

有关多项式的基本概念

对于求和式 $\sum a_n{x}^n$ ，如果是有限项相加，称为多项式，记作 $f(x)=\sum _{k=0}^n{a}_k{x}^k$。

一个多项式的度（也称次数）就是这个多项式最高次项的次数，记作 $\mathrm{deg}f$ 。

多项式的表示一般有两种方法：

1. 系数表示，就是形如 $f(x)=\sum _{k=0}^n{a}_k{x}^k$ 的形式。
2. 点值表示，即给出 $n+1$ 个点 $(x_0,y_0),(x_1,y_1),\dots ,(x_n,y_n)$ ，求一个 $n$ 次多项式使得这 $n+1$ 个点都在 $f(x)$ 上。

定义加法和乘法以后，多项式的基本运算和整数基本一致。（这个大家都会吧）

多项式的导数公式：

$${\left(\sum _{k=0}^{+\mathrm{\infty }}{a}_k{x}^k\right)}^{\prime }=\sum _{k=1}^{+\mathrm{\infty }}k{a}_k{x}^{k-1}$$

多项式的积分公式：

$${\int }_0^x\left(\sum _{k=0}^{+\mathrm{\infty }}{a}_k{x}^k\right)\text{d}x=\sum _{k=0}^{+\mathrm{\infty }}\frac{a_k{x}^{k+1}}{k+1}$$

一个长度为 $n+1$ 的序列 $a_i$ 的生成函数定义为 $a(x)=\sum _{k=0}^n{a}_k{x}^k$ 。（以下全文中 $f_i$ 的生成函数记为 $f(x)$ ）

多项式乘法

单位根（前置知识）

前置知识：复数

顾名思义，单位根就是满足 $x^n=1$ 的根。

众所周知，在复数域内， $n$ 次多项式有 $n$ 个根，因此 $x^n=1$ 也有 $n$ 个根，一般把第 $k$ 个根记作 ${\omega }_n^k$ 。

由于复数乘法的几何意义是模长相乘，辐角相加，因此 ${\omega }_n^k$ 的模长一定为 $1$ ，辐角为 $\frac{2\pi }{n}$ 的倍数。

为方便，我们设 ${\omega }_n^k$ 的辐角为 $\frac{2k\pi }{n}$ ，或者 $\frac{{360}^{\circ }k}{n}$ 。

单位根有三个重要的性质。对于任意正整数 $n$ 和整数 $k$：

$$\begin{array}{rl}{\omega }_n^n& =1\\ {\omega }_n^k& ={\omega }_{xn}^{xk}\\ {\omega }_{2n}^{k+n}& =-{\omega }_{2n}^k\end{array}$$

这三个性质在快速傅里叶变换中会被用到。

顺便一提，单位根还有一个比较重要的性质，是： $\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=0$ 。

DFT

又称离散傅里叶变换、Discrete Fourier Transform。

DFT的思想是将两个多项式的系数表示与点值表示互相转化。

因为多项式的系数表示直接相乘是 $O(n^2)$ 的，而点值表示相乘是 $O(n)$ 的。（只需要把每个点处的值相乘）

先考虑将系数表示转化为点值。

如果暴力去求，时间复杂度还是 $O(n^2)$ 的。（这个过程就被称为DFT）

那这个时候我们就要用到一个东西叫做：

FFT

又称快速傅里叶变换、Fast (Discrete) Fourier Transform。

我们先把多项式的次数 $n$ 增加至 $2^k-1$ 的形式。（一会儿我就会说明为什么要这样做了）

接下来，我们先尝试将系数表示转化为点值表示，这一步就被称之为FFT。

我们考虑一个 $n$ 次多项式有什么性质：

以 $n=7$ 为例：

$$f(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3+a_4{x}^4+a_5{x}^5+a_6{x}^6+a_7{x}^7$$

我们先按照次数的奇偶来分成两组，然后右边提出来一个 x：

$$\begin{array}{rl}f(x)& =(a_0+a_2{x}^2+a_4{x}^4+a_6{x}^6)+(a_{1}x+a_3{x}^3+a_5{x}^5+a_7{x}^7)\\ & =(a_0+a_2{x}^2+a_4{x}^4+a_6{x}^6)+x(a_1+a_3{x}^2+a_5{x}^4+a_7{x}^6)\end{array}$$

然后再用奇偶次次项数建立新的函数：

$$\begin{array}{rl}{f}_0(x)& =a_0+a_{2}x+a_4{x}^2+a_6{x}^3\\ f_1(x)& =a_1+a_{3}x+a_5{x}^2+a_7{x}^3\end{array}$$

那么原来的 $f(x)$ 就可以表示为：

$$f(x)=f_0\left(x^2\right)+x\times f_1\left(x^2\right)$$

于是我们可以发现一个性质：它可以分治！

也就是说，我们可以先将原来的式子的奇数项和偶数项分别处理（即分别计算FFT），求出原文中的 $f_0$ 和 $f_1$ ，再代入式子求解。

那怎么合并呢？

由于 $f_1$ 和 $f_2$ 都只有 $n/2$ 个点值，那么我们如果用以上那个式子的话，也只能推出 $f$ 的 $n/2$ 个点值。

但是，我们可以再考虑以下这个式子：

$$f(-x)=f_0\left(x^2\right)-x\times f_1\left(x^2\right)$$

这样的话，我们能用这个式子求出另外 $n/2$ 个点值。

你以为这就完了吗？还没有！

观察上面两个式子，我们发现：它只能求互为相反数的 $n/2$ 对 $x$ ，而如果在实数范围内考虑，右边就要求 $n/2$ 个正数的点值，显然行不通。

于是我们考虑单位根。

单位根有很好的的对称性，并且满足平方还是单位根这一特性。（因为 $({\omega }_n^k{)}^2={\omega }_{n/2}^k$ ）

如果我们把 $n$ 个单位根代入，那么式子就会变成这样：

$$\begin{array}{rl}f({\omega }_n^k)& =f_0({\omega }_{n/2}^k)+{\omega }_n^k{f}_1({\omega }_{n/2}^k)\\ f({\omega }_n^{k+n/2})& =f_0({\omega }_{n/2}^k)-{\omega }_n^k{f}_1({\omega }_{n/2}^k)\end{array}$$

（注意到之前我们让 $n$ 扩大到 $2^k-1$ 的形式，就是将项数扩大到 $2^k$ ，方便减半）

于是递归即可。

边界条件： $n=1$ 时，$f(x)=a_0$ 。

时间复杂度：$T(n)=2T(\frac{n}{2})+O(n)$ ，化简得 $O(n\log n)$ 。

逆FFT(IFFT)

现在我们已经解决了系数表示转换为点值表示，那么我们怎么将点值表示转换为系数表示呢？

我们知道有一个东西叫做高斯消元逆FFT。

它的证明需要依赖矩阵，这里不详细说明。（主要是作者也不是很会）

简单来说，你只需要把之前FFT的过程变成这样：

$$\begin{array}{rl}f({\omega }_n^k)& =f_0({\omega }_{n/2}^k)+{\omega }_n^{-k}{f}_1({\omega }_{n/2}^k)\\ f({\omega }_n^{k+n/2})& =f_0({\omega }_{n/2}^k)-{\omega }_n^{-k}{f}_1({\omega }_{n/2}^k)\end{array}$$

（就是系数变成它的共轭复数）

然后最终得出来的结果除以 $n$ 即可。

好的，到了这一步，我们已经会多项式乘法的基本操作了。

但是别着急，后面还有一个优化：

非递归写法

虽然 FFT 和 IFFT 都是 $O(n\log n)$ 的，但是它们有一个致命的缺陷：常数太大。

那我们应该怎么解决这个问题呢？

我们观察到 FFT 需要大量的复制数组的操作，这使得它的常数非常大。

有没有通过不移动数就能求出 FFT 的方法呢？

当然有。

我们先把整个过程压到一个数组中，再考虑之前的 FFT 都把数移到了哪些位置：（还是以 $n=7$ 举例，同一个数组的用中括号括起来）

下标	0	1	2	3	4	5	6	7
第一层	[0	1	2	3	4	5	6	7]
第二层	[0	2	4	6]	[1	3	5	7]
第三层	[0	4]	[2	6]	[1	5]	[3	7]
边界	[0]	[4]	[2]	[6]	[1]	[5]	[3]	[7]

可以看出最后一层每一位都是这个位置下标的二进制翻转过来。

而我们如果最开始把边界数组存下来，再将固定距离的数进行合并即可。

于是我们可以先预处理出每个数的二进制翻转后的结果 $re{v}_i$ ，顺便求出边界情况的数组 $b_i=a_{{rev}_i}$ 。

然后从最下面一层开始枚举，每层距离变为原来的 $2$ 倍，并对同一个区域内这个距离的数进行合并。（类似后缀排序的归并）

这样直接合并，时间复杂度仍为 $O(n\log n)$ ，但常数小了不少。

代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const double Pi=acos(-1);//用这个式子算pi要好一些
inline int read()
{
	char ch=getchar();int x=0,r=1;
	while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}
	while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return r?x:-x;
}
struct complexs
{
	double x,y;
	complexs(double xx=0,double yy=0){x=xx;y=yy;}
	complexs operator +(double b){return {x+b,y};}
	complexs operator -(double b){return {x-b,y};}
	complexs operator *(double b){return {x*b,y*b};}
	complexs operator /(double b){return {x/b,y/b};} 
	complexs operator +(complexs b){return {x+b.x,y+b.y};}
	complexs operator -(complexs b){return {x-b.x,y-b.y};}
	complexs operator *(complexs b){return {x*b.x-y*b.y,x*b.y+y*b.x};}
	complexs operator /(complexs b){return (complexs){x*b.x+y*b.y,-x*b.y+y*b.x}/(b.x*b.x+b.y*b.y);} 
};//虽然看上去很复杂但实际上就是简单的复数运算
typedef vector<complexs> vc;//vc就是vector<complexs>
int rev[2100000];//二进制翻转数组
void FFT(vc &x,int limit,bool type=1)//type=1是FFT，否则是IFFT
{
	for(int i=0;i<limit;++i)if(i<rev[i])swap(x[i],x[rev[i]]);//初始边界数组
	for(int len=1;len<limit;len<<=1)
	{
		complexs Wn(cos(Pi/len),(2*type-1)*sin(Pi/len));//单位根，注意IFFT的时候单位根要取反
		for(int i=0;i<limit;i+=2*len)//枚举每段左边界
		{
			complexs W(1,0),X,Y;
			for(int j=i;j<i+len;++j)
			{
				X=x[j];Y=W*x[j+len];
				x[j]=X+Y;x[j+len]=X-Y;//合并
				W=W*Wn;
			}
		}
	}
	if(!type)for(int i=0;i<=limit;++i)x[i]=x[i]/limit;//IFFT时每个数要除以n
}
int limit;
void operator *=(vc &a,vc b)
{
	int x=a.size()+b.size()-1;
	limit=1;
	while(1ull*limit<=a.size()-1+b.size()-1)limit<<=1;//注意limit是要大于等于项数，不是大于等于次数，项数=次数+1
	for(int i=0;i<limit;++i)rev[i]=(rev[i/2]/2+(i%2)*limit/2);//计算二进制翻转
	while(a.size()<1ull*limit)a.push_back(0);
	while(b.size()<1ull*limit)b.push_back(0);
	FFT(a,limit);FFT(b,limit);//系数转换成点值
	for(int i=0;i<limit;++i)a[i]=a[i]*b[i];//点值相乘
	FFT(a,limit,0);//点值转换成系数
	while(a.size()>1ull*x)a.pop_back();
}
vc a,b;
int n,m;
signed main()
{
	n=read();m=read();
	for(int i=0;i<=n;++i)a.push_back({1.0*read(),0}); 
	for(int i=0;i<=m;++i)b.push_back({1.0*read(),0});
	a*=b;
	for(complexs i:a)printf("%lld ",(int)(i.x+0.5));//精度丢失较大，建议四舍五入
	return 0;
}

三次变两次优化

省流：一般不用。

假设我们要求 $f(x)\ast g(x)$ ，考虑 $p(x)=f(x)+ig(x)$ 。

于是 $p^2(x)=f^2(x)-g^2(x)+2i\ast f(x)g(x)$ 。

发现它的虚部除以 $2$ 就是答案，输出即可。

NTT

又称快速数论变换、Number-Theoretic Transform。（在有的地方也被称为FNTT？总之都是这个意思）

虽然FFT好用，但它也有一些致命的弱点，例如复数乘法和三角函数常数太大，精度不高，不支持取模等等。

那我们能解决这些问题吗？不能

当然可以，于是我们就有了 NTT。

回忆一下FFT之所以用单位元是用了它的哪些性质。

• $\mathrm{\forall }0\le i,j<n,{\omega }_n^i\ne {\omega }_n^j$（确定点对个数时）
• ${\omega }_{kn}^{km}={\omega }_n^m$（倍增时）
• $({\omega }_n^m{)}^2=({\omega }_n^{m+n/2}{)}^2={\omega }_n^{2m}$（合并时）
• ${\omega }_1^0=1$（写起来方便）

于是我们发现数论中的原根也有这个性质。

前置知识：原根（你只要知道有这么个东西就行了，细节可以不用看）

更确切地说，我们设原根为 $g$ ，令 ${\omega }_n^m=g^{\frac{(p-1)m}{n}}$ 即可。

利用数论，它可以解决精度和取模问题。

其它代码基本和 FFT 一样。

注意：NTT不能用三次变两次优化！

证明

第一条由原根的定义，只要原根的指数不同那么值也不同。

由于 $\frac{p-1}{n}$ 总是整数（并且 $n$ 一定时是常数），且 $m$ 值不一致，故成立。

第二条：${\omega }_{kn}^{km}=g^{\frac{km(p-1)}{kn}}=g^{\frac{m(p-1)}{n}}={\omega }_n^m$

第三条：$({\omega }_n^m{)}^2=(g^{\frac{m(p-1)}{n}}{)}^2=g^{\frac{2m(p-1)}{n}}={\omega }_n^{2m}$

并且 $({\omega }_n^{m+n/2}{)}^2=(g^{\frac{(m+n/2)(p-1)}{n}}{)}^2=g^{\frac{2m(p-1)}{n}}\cdot g^{p-1}={\omega }_n^{2m}\cdot 1={\omega }_n^{2m}$

证毕。

代码（用vector实现）

#include<bits/stdc++.h>
using namespace std;
#define int long long
typedef vector<int> vi;//vi指vector<int> 
const int mod=998244353,g=3,ig=332748118;//ig是g的逆元
inline int read()
{
	char ch=getchar();int x=0,r=1;
	while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}
	while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return r?x:-x;
}
int mul(int x,int b)//快速幂
{
	int ans=1;
	while(b)
	{
		if(b&1)ans=(ans*x)%mod;
		x=(x*x)%mod;
		b>>=1;
	}
	return ans;
}
int rev[2100000];
void NTT(vi &x,int limit,bool type=1)
{
	for(int i=0;i<limit;++i)if(i<rev[i])swap(x[i],x[rev[i]]);
	for(int len=1;len<limit;len<<=1)
	{
		int Wn=mul(type?g:ig,(mod-1)/(2*len));//利用原根代替FFT中的单位元
		for(int i=0;i<limit;i+=2*len)
		{
			int W=1,X,Y;
			for(int j=i;j<i+len;++j)
			{
				X=x[j];Y=W*x[j+len]%mod;
				x[j]=(X+Y)%mod;x[j+len]=(X-Y+mod)%mod;//合并同FFT，但是要加上取模运算
				W=(W*Wn)%mod;
			}
		}
	}
	if(!type)
	{
		int invs=mul(limit,mod-2);
		for(int i=0;i<limit;++i)x[i]=(x[i]*invs)%mod;
	}
}
int limit;
void operator *=(vi &a,vi b)
{
	int x=a.size()+b.size()-1;
	limit=1;
	while(1ull*limit<=a.size()-1+b.size()-1)limit<<=1;
	for(int i=0;i<limit;++i)rev[i]=(rev[i/2]/2+(i%2)*limit/2);
	while(a.size()<1ull*limit)a.push_back(0);
	while(b.size()<1ull*limit)b.push_back(0);
	NTT(a,limit);NTT(b,limit);//FFT变为NTT
	for(int i=0;i<limit;++i)a[i]=(a[i]*b[i])%mod;
	NTT(a,limit,0);
	while(a.size()>1ull*x)a.pop_back();
}
vi a,b;
int n,m;
signed main()
{
	n=read();m=read();
	for(int i=0;i<=n;++i)a.push_back(read());
	for(int i=0;i<=m;++i)b.push_back(read());
	a*=b;
	for(int i:a)printf("%lld ",i);
	return 0;
}

MTT

省流：解决任意模数的问题，一般不会用，但出现时比较恶心，可以写多次卷积+CRT（中国剩余定理），就是常数容易爆。

我们发现NTT只适用于模数为 $a\ast 2^n+1$ 的形式。（2^n要大于多项式项数，否则也会炸）

那如果模数不是 $998244353$ 这样的和谐的数字，是 ${10}^9+7$ 之类的呢？

解法一

假设原来的多项式值域是 $m={10}^9$ ，个数是 $n={10}^5$ （一般都是这个范围），那么能产生的最大数也不过就是 $m{n}^2={10}^{26}$ 。

考虑对多个模数卷积（ $3$ 个就够了），再用CRT（中国剩余定理）合并， $9$ 次FFT，能过一般题，但容易炸。

解法二

考虑拆系数。

我们先顺手把每个系数先模上 $mod$ ，于是系数的值域一定 $<mod$ 。

把系数拆成 $c=a\ast 2^{15}+b$ 的形式。

把 $a$ 和 $b$ 拆成两个多项式，相乘的值域为 $n\ast \sqrt{m}\ast \sqrt{m}={10}^{14}$

于是 $c_1\ast c_2=(a_1\ast 2^{15}+b_1)(a_2\ast 2^{15}+b_2)=a_1{a}_2\ast 2^{30}+(a_1{b}_2+a_2{b}_1)\ast 2^{15}+b_1{b}_2$ 。

对这个式子的 $a_1,a_2,b_1,b_2$ 分别做一次FFT，$a_1{a}_2,a_1{b}_2,a_2{b}_1,b_1{b}_2$ 分别做一次IFFT，总共要用 $8$ 次，常数还是比较大。

并且我们不能用可爱的NTT了，这意味着我们要面对FFT的精度问题（毕竟上界是 ${10}^{1}4$ ），建议写 $\text{long double}$ 。

以及，千万不要相信自己背 $\pi$ 的能力，建议写 $\arccos (-1)$ 。

解法三

我们考虑在解法二的基础上进行优化。

我们现在有四个多项式 $a_1,a_2,b_1,b_2$ ，我们要求 $a_1{a}_2,a_1{b}_2,a_2{b}_1,b_1{b}_2$ 的值。

上面的三次变两次优化告诉我们：有的时候塞一个复多项式进去可以把常数变得更小。

不妨设 $p=a_1+i{a}_2$ ， $q=b_1+i{b}_2$ 。

则 $pq=a_1{b}_1-a_2{b}_2+i(a_1{b}_2+a_2{b}_1)$ 。

再取 $p^{\prime }=a_1-i{a}_2$ ， 则 $p^{\prime }q=a_1{b}_1+a_2{b}_2+i(a_1{b}_2-a_2{b}_1)$。

因此，我们对 $p,p^{\prime },q$ 做FFT，对 $pq,p^{\prime }q$ 做IFFT，剩下的就是一些基本的多项式加减操作了，总共用 $5$ 次。

还是要写 $\text{long double}$ 。

多项式的各种运算

多项式乘法逆（倍增法）

多项式乘法逆的定义：对于一个 $n$ 次多项式 $f(x)$ 来说， $f(x)$ 的逆是满足 $f(x)\cdot g(x)\equiv 1(\mathrm{mod}{x}^{n+1})$ 的多项式 $g(x)$ 。

我们还是考虑递归。

首先便于计算，我们还是把 $n$ 扩大到 $2^k-1$ 的形式。（有些题解里不是这样写的，你也可以直接递归，但这样方便一些）

由于高次项的系数不会影响低次项的逆元，我们这步“扩大”可以直接补 $0$ 。

为了下文叙述方便，我们记 $m=x^{n+1},m^{\prime }=x^{\frac{n+1}{2}}$ 。

我们可以先求出 $f(x)$ 对 $m^{\prime }$ 的逆，记作 $f_0^{-1}(x)$ 。

于是我们有：

$$\begin{array}{rlr}f(x)f_0^{-1}(x)& \equiv 1& (\mathrm{mod}{m}^{\prime })\\ f(x)f^{-1}(x)& \equiv 1& (\mathrm{mod}{m}^{\prime })\\ f^{-1}(x)-f_0^{-1}(x)& \equiv 0& (\mathrm{mod}{m}^{\prime })\end{array}$$

将最后一个式子平方可得：

$$f^{-2}(x)-2f^{-1}(x)f_0^{-1}(x)+f_0^{-2}(x)\equiv 0(\mathrm{mod}m)$$

两边同乘 $f(x)$ ：

$$f^{-1}(x)-2f_0^{-1}(x)+f(x)f_0^{-2}(x)\equiv 0(\mathrm{mod}m)$$

移项得：

$$f^{-1}\left(x\right)\equiv f_0^{-1}\left(x\right)(2-f\left(x\right)f_0^{-1}\left(x\right))(\mathrm{mod}m)$$

倍增即可。

时间复杂度： $T(n)=T(\frac{n}{2})+O(n\log n)$ ，即 $T(n)=O(n\log n)$ 。

代码

//没错我又更了一堆插件
void read(int n,vi &a){for(int i=0;i<=n;++i)a.push_back(read());}//输入，n为次数，或者项数+1
void print(vi &a){for(int i:a)printf("%lld ",i);puts("");}//输出
void cut(vi &a,int x){while(a.size()>1ull*x)a.pop_back();while(a.size()<1ull*x)a.push_back(0);}//将a的系数个数强制设为x，注意系数个数和deg是两个不同的东西
void copy(vi &a,int x,vi &b)//将a的前x个系数复制到b，不足补0
{
	if(1ull*x<=a.size()){b.clear();for(int i=0;i<x;++i)b.push_back(a[i]);}
	else {b=a;while(b.size()<1ull*x)b.push_back(0);}
}
void operator +=(vi &a,int x){a.front()=(a.front()+x)%mod;}//多项式+常数
void operator +=(vi &a,vi &b)//多项式相加
{
	while(a.size()<b.size())a.push_back(0);
	for(int i=0;1ull*i<b.size();++i)a[i]=(a[i]+b[i])%mod;
}
void operator --(vi &a){for(int i=0;1ull*i<a.size();++i)a[i]=(mod-a[i])%mod;}//将多项式取负，相当于0-a 
void operator -=(vi &a,int x){a.front()=(a.front()-x+mod)%mod;}//多项式-常数
void operator -=(vi &a,vi &b)//多项式相减
{
	while(a.size()<b.size())a.push_back(0);
	for(int i=0;1ull*i<b.size();++i)a[i]=(a[i]-b[i]+mod)%mod;
}
vi inv(vi &a)//逆多项式部分
{
	int x=1;
	vi ans={mul(a.front(),mod-2)},b,now;//显然a对于x^1的逆就是a0的逆元
	while(1ull*x<a.size())
	{
		x<<=1;//x为系数个数
		copy(a,x,b);//将a的后几位存下来，避免ans每次直接与a相乘，影响时间复杂度
		now=ans;ans*=b;ans-=2;--ans;ans*=now;//按式子进行计算，注意代码中--的定义与整数不同
		cut(ans,x);//去除多余位
	}
	cut(ans,a.size());//让ans次数和a一样
	return ans;
}

多项式除法

顾名思义，就是一个多项式 $F(x)$ 除以另一个多项式 $G(x)$ 的答案。

我们直接对 $G(x)$ 取逆，再乘上 $F(x)$ 即可。

好吧，真正的多项式除法实际上是带鱼带余除法。

或者说得正式点：给定一个 $n$ 次多项式 $f(x)$ 和一个 $m$ 次多项式 $g(x)$ ，请求出多项式 $q(x)$, $r(x)$，满足以下条件：

• $\mathrm{deg}q(x)=n-m$，$\mathrm{deg}r(x)=m-1$
• $f(x)=q(x)\ast g(x)+r(x)$

要解决这样一个问题，我们首先构造这样一个变换：

$$f^R(x)=x^{\mathrm{deg}f}f(\frac{1}{x})$$

观察到这个变换的实质是反转 $f(x)$ 的系数，即 $\text{vector}$ 当中的 $\text{reverse}$ 函数。

并且实际上有 $(f^R{)}^R(x)=f(x)$ ，因此我们只需要求出 $f^R(x)$ 就能求出 $f(x)$ 。

（注： $f^R(x)$ 和 $f^{-1}(x)$ 是两个不同的东西，一个是反转系数，一个是乘法逆）

由 $f(x)=q(x)\ast g(x)+r(x)$ 可知：

$$f(\frac{1}{x})=q(\frac{1}{x})\ast g(\frac{1}{x})+r(\frac{1}{x})$$

两边同时乘上 $x^n$ 得：

$$f^R(x)=q^R(x)\ast g^R(x)+x^{n-m+1}{r}^R(x)$$

注意到我们只需要求出 $q(x)$ 和 $r(x)$ 中的任意一个就能求出另一个。

我们发现上式中 $\mathrm{deg}q(x)=n-m$ ，而 $r^R(x)$ 的系数为 $x^{n-m+1}$ 。

因此我们把式子两边模上 $x^{n-m+1}$ :

$$f^R(x)\equiv g^R(x)q^R(x)(\mathrm{mod}{x}^{n-m+1})$$

$$f^R(x)(g^R(x){)}^{-1}\equiv q^R(x)(\mathrm{mod}{x}^{n-m+1})$$

于是多项式求逆解决 $q(x)$ ，再代入原式解出 $r(x)$ 即可。

注意要先把 $g^R(x)$ 的次数设为 $n-m$ 再求逆。

时间复杂度 $O(n\log n)$ 。

代码

typedef pair<vi,vi> pvv;//除法返回答案时要用pair，你也可以设几个全局变量
#define mk make_pair
#define F first
#define S second
void R(vi &a){reverse(a.begin(),a.end());}//简单的reverse
pvv operator /(vi a,vi b)
{
	vi q,r;
	q=b;R(q);cut(q,a.size()-b.size()+1);//把q设为bR，把q的次数设为n-m（size设为n-m+1）
	r=a;R(r);//r作为临时变量存储aR
	q=inv(q);q*=r;//计算qR
	cut(q,a.size()-b.size()+1);R(q);//将qR多余位数截掉后还原q
	r=q;r*=b;r-=a;--r;cut(r,b.size()-1);//计算r
	return mk(q,r);
}

多项式开方

给定一个 $n$ 次多项式 $f(x)$ ，请求出多项式 $g(x)$ ，满足 $g^2(x)=f(x)$ 。

设 $f(x)$ 的系数为 $a_0,a_1,\cdots ,a_n$ ， $g(x)$ 的系数为 $b_0,b_1,\cdots ,b_n$ 。

首先我们先求出 $g(x)$ 的常数项 $b_0$ 。

代入 $x=0$ 可得： $b_0^2=a_0$ ，因此如果 $a_0$ 在模意义下没有平方根就无解。

于是我们取 $b_0=\sqrt{a_0}$ 。

（由于 $a_0$ 的平方根可能不止一个，我们随便取一个值作为 $b_0$ ；选择不同的 $b_0$ 会求出不同的 $b(x)$ 。）

以下我们先假定 $a_0=1$ ，因为这样讨论起来比较方便。

老规矩，我们继续把 $n$ 扩大到 $2^k-1$ 的形式。

假设我们求出了 $f(x)$ 在模 $x^{\frac{n}{2}}$ 意义下的平方根 $g_0(x)$ 。

记 $m=x^n,m^{\prime }=x^{\frac{n}{2}}$ 。则有：

$$f(x)-g_0^2(x)\equiv 0(\mathrm{mod}{m}^{\prime })$$

两边同时平方，得：

$$\begin{array}{rl}(f(x)-g_0^2(x){)}^2& \equiv 0(\mathrm{mod}m)\\ (f(x)+g_0^2(x){)}^2& \equiv 4f(x)g_0^2(x)(\mathrm{mod}m)\end{array}$$

将 $4g_0^2(x)$ 移到左边得：

$${\left(\frac{f(x)+g_0^2(x)}{2g_0(x)}\right)}^2\equiv f(x)(\mathrm{mod}m)$$

两边同时开方：

$$\begin{array}{rl}\frac{f(x)+g_0^2(x)}{2g_0(x)}& \equiv g(x)(\mathrm{mod}m)\\ 2^{-1}(f(x)g_0^{-1}(x)+g_0(x))& \equiv g(x)(\mathrm{mod}m)\end{array}$$

倍增计算即可。

时间复杂度： $T(n)=T(\frac{n}{2})+O(n\log n)$ ，即 $T(n)=O(n\log n)$ 。

代码

vi sqrt(vi &a)
{
	int x=1;
	vi ans={1},b,now;
	while(1ull*x<a.size())
	{
		x<<=1;
		cut(ans,x);//这里求逆的时候要
		copy(a,x,b);
		now=ans;ans=inv(ans);ans*=b;ans+=now;ans*=i2;
		cut(ans,x);
	}
	cut(ans,a.size());
	return ans;
}

（为了直观，我写的常数比较大，建议读者写的时候不要每次都做一遍cut和copy，建议多传几个参数）

然后我们解决 $a_0$ 为任意数时的问题。（一般不会这么毒瘤）

若 $a_0\ne 0$ ，我们用二次剩余即可。

否则，我们考虑原式中因式 $x$ 的次数 $k$ 。（即：找到一个最大的 $k$ ，使得 $a_0,a_1,a_2,\cdots ,a_k=0$ ）

如果 $k$ 是奇数，则无解；如果 $k$ 是偶数，输出 $\sqrt{\frac{a}{x^k}}\ast x^{\frac{k}{2}}$ 即可。

分治FFT

经典例题：给你一个数组 $g_i$ ，求数组 $f_i$ ，使得 $f_i=\sum _{j=1}^i{f}_{i-j}{g}_j$ 且 $f_0=1$。

我们先来对比一下普通卷积和分治FFT的式子：

普通卷积：$h_i=\sum _{j=1}^i{f}_{i-j}{g}_j$

分治FFT：$f_i=\sum _{j=1}^i{f}_{i-j}{g}_j$

我们发现：分治FFT实际上和普通卷积比较像，就是输出的数组就是它本身。

但这个性质好像没啥卵用

我们发现，对于每个区间 $[l,r]$ ，记 $mid=\lfloor \frac{l+r}{2}\rfloor$ ，考虑一种类似于CDQ分治的算法：

1. 递归 $[l,mid]$ 算出递归区间内的答案。
2. 把 $[l,mid]$ 的贡献加到 $[mid+1,r]$ 上。
3. 递归 $[mid+1,r]$ 。

其中第二步我们再详细讲一下。

假设当前区间为 $[0,5]$ ，则我们需要

$$\begin{array}{rl}{f}_3& +=f_2\ast g_1+f_1\ast g_2+f_0\ast g_3\\ f_4& +=f_2\ast g_2+f_1\ast g_3+f_0\ast g_4\\ f_5& +=f_2\ast g_3+f_1\ast g_4+f_0\ast g_5\end{array}$$

于是，$f_n+=\sum _{i=l}^{mid}{f}_i{g}_{n-i}=[x^n]f^{\prime }(x)g^{\prime }(x)$ ，其中 $f^{\prime }(x)$ 代表的是只考虑 $[l,mid]$ 这一段时的 $f(x)$ ， $g^{\prime }(x)$ 代表只考虑 $[1,r-l]$ 时的 $g(x)$ 。

卷积即可。

时间复杂度： $T(n)=2T(\frac{n}{2})+O(n\log n)$ ，即 $T(n)=O(n{\log }^{2}n)$ 。

代码

void copy(vi &a,int l,int r,vi &b)//把a的第l位到第r位复制给b
{
	b.clear();
	for(int i=l;i<=r;++i)b.push_back(a[i]);
}
void div_FFT(vi &a,vi &b,int l=-1,int r=-1)
{
	if(l==-1&&r==-1)
	{
		cut(a,b.size());
		l=0;r=b.size()-1;//初始化
	}
	if(l==r)return;
	int mid=(l+r)>>1;
	div_FFT(a,b,l,mid);//分治左区间
	vi c,d;
	copy(b,1,r-l,c);copy(a,l,mid,d);
	c*=d;
	for(int i=mid+1;i<=r;++i)a[i]=(a[i]+c[i-l-1])%mod;//卷积处理对右半的贡献
	div_FFT(a,b,mid+1,r);//分治右区间
}

但是先别走！这道题也可以用多项式求逆解决！！！

我们发现

$$f(x)g(x)\equiv \sum _{i=0}^{n-1}{x}^i\sum _{j=0}^i{f}_{i-j}{g}_j(\mathrm{mod}{x}^n)$$

由于 $g_0=0$，所以

$$f(x)g(x)\equiv \sum _{i=1}^{n-1}{x}^i\sum _{j=1}^i{f}_{i-j}{g}_j(\mathrm{mod}{x}^n)$$

并且

$$\begin{array}{rl}f(x)& \equiv \sum _{i=0}^{n-1}{x}^i{f}_i\\ & \equiv x^0{f}_0+\sum _{i=1}^{n-1}{x}^i{f}_i\\ & \equiv f_0+\sum _{i=0}^{n-1}{x}^i\sum _{j=1}^i{f}_{i-j}{g}_j(\mathrm{mod}{x}^n)\end{array}$$

则有 $f(x)g(x)\equiv f(x)-f_0$ ，于是就有 $f(x)(1-g(x))\equiv f_0$，从而 $f(x)\equiv f_0(1-g(x){)}^{-1}$，多项式求逆即可。

时间复杂度 $O(n\log n)$。

多项式 $\ln$

各位都知道整数的 $\ln$ 和 $\exp$ 吧。（ $exp(x)$ 实际上就是 $e^x$）

后面的证明需要依赖微积分，不会的同学请自行跳过。

实际上， $\ln$ 和 $\exp$ 的完整定义在这：

$$\ln f(x)=-\sum _{i=1}^{\mathrm{\infty }}\frac{(1-f(x){)}^i}{i}$$

$$\exp f(x)=\sum _{i=0}^{\mathrm{\infty }}\frac{f(x{)}^i}{i!}$$

但是，直接用定义式时间复杂度会炸。

首先，对于多项式 $f(x)$，若 $\ln (x)$ 存在，则有 $[x_0]f(x)=1$ 。（否则常数项不收敛）

我们知道 $\ln x$ 的导数为 $\frac{1}{x}$ ，不妨设 $\ln f(x)=g(x)$ ，两边同时对 $x$ 求导可得：

$$\frac{1}{f(x)}\cdot f^{\prime }(x)=g^{\prime }(x)$$

$$\frac{f^{\prime }(x)}{f(x)}=g^{\prime }(x)$$

再同时积分： $$g(x)=\int \left(\frac{f^{\prime }(x)}{f(x)}\right)$$

因此我们只需要多项式求导、积分、求逆和乘积就可以了。

时间复杂度 $O(n\log n)$。

代码

int invv[2100000];
void init(int n=2099999)
{
	invv[1]=1;
	for(int i=2;i<=n;++i)invv[i]=(mod-mod/i)*invv[mod%i]%mod;
}//线性求逆元
void D(vi &a){for(int i=1;1ull*i<a.size();++i)a[i-1]=i*a[i]%mod;a.pop_back();}//求导
void I(vi &a){a.push_back(0);for(int i=a.size()-1;i>=1;--i)a[i]=invv[i]*a[i-1]%mod;a[0]=0;}//积分
void ln(vi &a)
{
	vi b=inv(a);
	D(a);
	a*=b;I(a);
	cut(a,b.size());
}

多项式 $\exp$

首先，对于多项式 $f(x)$，若 $\exp (x)$ 存在，则有 $[x_0]f(x)=0$ 。（否则常数项不收敛）

我们知道 $\exp x$ 的导数就是它本身，不妨设 $\exp f(x)=g(x)$ 。

两边求导得： $g(x)\ast f^{\prime }(x)=g^{\prime }(x)$

提取系数后就变成了 $g_n^{\prime }=\sum _{i=0}^n{g}_{n-i}{f}_i^{\prime }$

由多项式的求导，有 $(n+1)g_{n+1}=\sum _{i=0}^n{g}_{n-i}{f}_{i+1}(i+1)$

即 $n{g}_n=\sum _{i=1}^n(i\cdot g_{n-i}{f}_i)$，$g_n=\frac{1}{n}\sum _{i=1}^n(i\cdot g_{n-i}{f}_i)$。

并且 $g_0=1$ ，分治FFT即可。

虽然和分治FFT的模板不同，但实际上差不多，我们把 $i$ 乘到 $f_i$ 上，再在 $l=r$ 处把常数 $\frac{1}{n}$ 乘上。

其实，把上述式子再变换也可以求出 $\ln$ 的递推形式，这里不再阐述，感兴趣的同学可以自己推一下。

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

代码

void div_FFT(vi &a,vi &b,bool isexp=0,int l=-1,int r=-1)//对之前的分治FFT做了点手脚
{
	if(l==-1&&r==-1)
	{
		if(a.empty())a.push_back(1);
		cut(a,b.size());
		l=0;r=b.size()-1;
	}
	if(l==r){if(isexp&&l)a[l]=(a[l]*invv[l])%mod;return;}//边界条件略改一下
	int mid=(l+r)>>1;
	div_FFT(a,b,isexp,l,mid);
	vi c,d;
	copy(b,1,r-l,c);copy(a,l,mid,d);
	c*=d;
	for(int i=mid+1;i<=r;++i)a[i]=(a[i]+c[i-l-1])%mod;
	div_FFT(a,b,isexp,mid+1,r);
}
vi exp(vi &a)
{
	for(int i=0;1ull*i<a.size();++i)a[i]=(a[i]*i)%mod;
	vi b={};
	div_FFT(b,a,1);//1代表分治FFT最终的系数要除以下标
	return b;
}

多项式快速幂

我们考虑封装多项式乘法，再倍增即可。

好吧，正解是这样的：

我们先保证 $a_0=1$ 。

众所周知，多项式的 $\ln$ 有个性质： $\ln f^k(x)=k\ln f(x)$ 。

于是我们输出 $\exp (k\cdot \ln f(x))$ 即可。

时间复杂度 $O(n\log n)$ 。

但这还没完！

我们发现洛谷上的模板题是这样的：

$0<k\le {10}^{{10}^5}$

那怎么办呢？还要写个高精？

肯定不可能，我们需要知道一个性质：

$$f(x{)}^k\equiv f(x{)}^{kmodp}(\mathrm{mod}{x}^n)$$

这个式子在 $n<p$ 时成立，而 $n\ge p$ 的情况几乎不考，所以就做完了。

证明

我们先证明一下这样一个结论：$f(x{)}^p\equiv f(x^p)$。

显然有个式子叫做 $(a+b{)}^p\equiv a^p+b^p$ 。（组合数拆分系数）

不妨设上述结论对所有 $k-1$ 次多项式成立。

设 $f(x)$ 的次数为 $k$ ，并且 $f(x)=a_k{x}^k+g(x)$，

则 $f(x{)}^p\equiv g(x{)}^p+(a_k{x}^k{)}^p\equiv g(x^p)+a_k^p{x}^{pk}\equiv g(x^p)+a_k{x}^{pk}=f(x^p)$。

有了这个式子以后，由于 $n<p$ ，所以 $f(x{)}^{p}mod{x}^n\equiv f(x^p)mod{x}^n\equiv a_0{x}^0\equiv 1$ ，证毕。

有些同学可能到了这里还有一个疑问：为什么整数快速幂的时候指数要模 $p-1$ ，但多项式快速幂就变成模 $p$ 了呢？

换句话说，如果把 $f(x)$ 带入具体的值，那么模数应该模上 $p-1$ ，而不是 $p$ 。

其实，这两个式子还是有一点差别的。

快速幂的式子是：

$$f(x{)}^k\equiv f(x{)}^{kmod(p-1)}(\mathrm{mod}p)$$

而在 $n<p$ 的条件下，多项式快速幂的式子是：

$$f(x{)}^{k}mod{x}^n\equiv f(x{)}^{kmodp}mod{x}^n$$

也就是说一个式子算完了，而另一个式子只计算了前 $n$ 项，因此它们推导的方式也不同。

一定要注意多项式快速幂的证明要依赖 $n<p$ 这个限制。

有了这个性质，我们就能很轻松地解决刚开始的问题了，只需要在读入的时候将 $k$ 模上 $p$ 即可，根本不需要高精。

啥？怎么在读入的时候将 $k$ 模上 $p$ ？你想想快读是怎么读入的就知道了。

啥？还不会？快读每次乘 $10$ 都取模啊！

代码

inline int read(int mod=0)//修改后的快读，加入了取模
{
	char ch=getchar();int x=0,r=1;
	while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';if(mod)x%=mod;ch=getchar();}//这里有改动
	return r?x:-x;
}
void mul(vi &a,int k)
{
	ln(a);
	a*=k;
	a=exp(a);
}

接下来我们考虑 $a_0\ne 1$ 的情况。

当 $a_0\ne 0$ 时，我们先把整个式子除以 $a_0$ ，然后再进行快速幂，最后再把整个式子乘上 $a_0^k$ 即可。

对于 $a_0=0$ 的情况，我们就把原式中的因式 $x^m$ 提出来，转化成上一种情况求解。

注意这里要同时保留 $kmodp$ 和 $kmod(p-1)$ 的值，因此不能直接快读，建议把结果存到字符串里。

多项式全家桶代码

由于之前的代码段太多了，函数也比较多，所以这里整合了一下。

建议大家先把所有的函数写完再一起卡常，防止自己看不懂

全家桶代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
typedef vector<int> vi;
typedef pair<vi,vi> pvv;
#define mk make_pair
#define F first
#define S second
const int mod=998244353,g=3,ig=332748118,i2=499122177;
inline int read(int mod=0)
{
	char ch=getchar();int x=0,r=1;
	while(ch<'0'||ch>'9'){if(ch=='-')r=0;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+ch-'0';if(mod)x%=mod;ch=getchar();}
	return r?x:-x;
}
int mul(int x,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1)ans=(ans*x)%mod;
		x=(x*x)%mod;
		b>>=1;
	}
	return ans;
}
int invv[2100000];
void init(int n=2099999)
{
	invv[1]=1;
	for(int i=2;i<=n;++i)invv[i]=(mod-mod/i)*invv[mod%i]%mod;
}
int rev[2100000];
void read(int n,vi &a,int l=0){for(int i=0;i<=l+n-1;++i)a.push_back(i<l?0:read());}
void print(vi &a){for(int i:a)printf("%lld ",i);puts("");}
void D(vi &a){for(int i=1;1ull*i<a.size();++i)a[i-1]=i*a[i]%mod;a.pop_back();}
void I(vi &a){a.push_back(0);for(int i=a.size()-1;i>=1;--i)a[i]=invv[i]*a[i-1]%mod;a[0]=0;}
void NTT(vi &x,int limit,bool type=1)
{
	for(int i=0;i<limit;++i)if(i<rev[i])swap(x[i],x[rev[i]]);
	for(int len=1;len<limit;len<<=1)
	{
		int Wn=mul(type?g:ig,(mod-1)/(2*len));
		for(int i=0;i<limit;i+=2*len)
		{
			int W=1,X,Y;
			for(int j=i;j<i+len;++j)
			{
				X=x[j];Y=W*x[j+len]%mod;
				x[j]=(X+Y)%mod;x[j+len]=(X-Y+mod)%mod;
				W=(W*Wn)%mod;
			}
		}
	}
	if(!type)
	{
		int invs=mul(limit,mod-2);
		for(int i=0;i<limit;++i)x[i]=(x[i]*invs)%mod;
	}
}
int limit;
void cut(vi &a,int x){while(a.size()>1ull*x)a.pop_back();while(a.size()<1ull*x)a.push_back(0);}
void copy(vi &a,int x,vi &b)
{
	if(1ull*x<=a.size()){b.clear();for(int i=0;i<x;++i)b.push_back(a[i]);}
	else {b=a;while(b.size()<1ull*x)b.push_back(0);}
}
void copy(vi &a,int l,int r,vi &b)
{
	b.clear();
	for(int i=l;i<=r;++i)b.push_back(a[i]);
}
void operator +=(vi &a,int x){a.front()=(a.front()+x)%mod;}
void operator +=(vi &a,vi &b)
{
	while(a.size()<b.size())a.push_back(0);
	for(int i=0;1ull*i<b.size();++i)a[i]=(a[i]+b[i])%mod;
}
void operator -=(vi &a,int x){a.front()=(a.front()-x+mod)%mod;}
void operator -=(vi &a,vi &b)
{
	while(a.size()<b.size())a.push_back(0);
	for(int i=0;1ull*i<b.size();++i)a[i]=(a[i]-b[i]+mod)%mod;
}
void operator *=(vi &a,int x){for(int i=0;1ull*i<a.size();++i)a[i]=((a[i]*x)%mod+mod)%mod;}
void operator *=(vi &a,vi b)
{
	int x=a.size()+b.size()-1;
	limit=1;
	while(1ull*limit<=a.size()-1+b.size()-1)limit<<=1;
	for(int i=0;i<limit;++i)rev[i]=(rev[i/2]/2+(i%2)*limit/2);
	while(a.size()<1ull*limit)a.push_back(0);
	while(b.size()<1ull*limit)b.push_back(0);
	NTT(a,limit);NTT(b,limit);
	for(int i=0;i<limit;++i)a[i]=(a[i]*b[i])%mod;
	NTT(a,limit,0);
	cut(a,x);
}
vi inv(vi &a)
{
	int x=1;
	vi ans={mul(a.front(),mod-2)},b,now;
	while(1ull*x<a.size())
	{
		x<<=1;
		copy(a,x,b);
		now=ans;ans*=b;ans*=-1;ans+=2;ans*=now;
		cut(ans,x);
	}
	cut(ans,a.size());
	return ans;
}
void R(vi &a){reverse(a.begin(),a.end());}
pvv operator /(vi a,vi b)
{
	vi q,r;
	q=b;R(q);cut(q,a.size()-b.size()+1);
	r=a;R(r);
	q=inv(q);q*=r;
	cut(q,a.size()-b.size()+1);R(q);
	r=q;r*=b;r-=a;r*=-1;cut(r,b.size()-1);
	return mk(q,r);
}
vi sqrt(vi &a)
{
	int x=1;
	vi ans={1},b,now;
	while(1ull*x<a.size())
	{
		x<<=1;
		cut(ans,x);copy(a,x,b);
		now=ans;ans=inv(ans);ans*=b;ans+=now;ans*=i2;
		cut(ans,x);
	}
	cut(ans,a.size());
	return ans;
}
void div_FFT(vi &a,vi &b,bool isexp=0,int l=-1,int r=-1)
{
	if(l==-1&&r==-1)
	{
		if(a.empty())a.push_back(1);
		cut(a,b.size());
		l=0;r=b.size()-1;
	}
	if(l==r){if(isexp&&l)a[l]=(a[l]*invv[l])%mod;return;}
	int mid=(l+r)>>1;
	div_FFT(a,b,isexp,l,mid);
	vi c,d;
	copy(b,1,r-l,c);copy(a,l,mid,d);
	c*=d;
	for(int i=mid+1;i<=r;++i)a[i]=(a[i]+c[i-l-1])%mod;
	div_FFT(a,b,isexp,mid+1,r);
}
void ln(vi &a)
{
	vi b=inv(a);
	D(a);
	a*=b;I(a);
	cut(a,b.size());
}
vi exp(vi &a)
{
	for(int i=0;1ull*i<a.size();++i)a[i]=(a[i]*i)%mod;
	vi b={};
	div_FFT(b,a,1);
	return b;
}
void mul(vi &a,int k)
{
	ln(a);
	a*=k;
	a=exp(a);
}
signed main()
{
	
	return 0;
}

多项式牛顿迭代

开个坑，以后再写。