多项式：从入门到全家桶

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

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

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

有关多项式的基本概念

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

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

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

系数表示，就是形如 $f(x)=\sum_{k=0}^n a_kx^k$ 的形式。
点值表示，即给出 $n+1$ 个点 $(x_{0}, y_{0}), (x_{1}, y_{1}), \dots, (x_{n}, y_{n})$ ，求一个 $n$ 次多项式使得这 $n+1$ 个点都在 $f(x)$ 上。

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

多项式的导数公式：

\[\left(\sum_{k=0}^{+\infty}a_kx^k\right)'=\sum_{k=1}^{+\infty}ka_kx^{k-1} \]

多项式的积分公式：

\[\int_0^x\left(\sum_{k=0}^{+\infty}a_kx^k\right)\text{d}x=\sum_{k=0}^{+\infty}\frac{a_kx^{k+1}}{k+1} \]

一个长度为 $n+1$ 的序列 $a_i$ 的生成函数定义为 $a(x)=\sum_{k=0}^n a_kx^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{aligned} \omega_n^n&=1\\ \omega_n^k&=\omega_{xn}^{xk}\\ \omega_{2n}^{k+n}&=-\omega_{2n}^k\\ \end{aligned} \]

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

顺便一提，单位根还有一个比较重要的性质，是： $\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_1x + a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 \]

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

\[\begin{aligned} f(x) &= (a_0+a_2x^2+a_4x^4+a_6x^6) + (a_1x+a_3x^3+a_5x^5+a_7x^7)\\ &= (a_0+a_2x^2+a_4x^4+a_6x^6) + x(a_1+a_3x^2+a_5x^4+a_7x^6) \end{aligned} \]

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

\[\begin{aligned} f_0(x) &= a_0+a_2x+a_4x^2+a_6x^3\\ f_1(x) &= a_1+a_3x+a_5x^2+a_7x^3 \end{aligned} \]

那么原来的 $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{aligned} 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{aligned} \]

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

于是递归即可。

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

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

逆FFT(IFFT)

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

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

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

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

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

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

然后最终得出来的结果除以 $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]

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

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

于是我们可以先预处理出每个数的二进制翻转后的结果 $rev_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)*g(x)$ ，考虑 $p(x)=f(x)+ig(x)$ 。

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

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

NTT

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

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

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

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

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

$\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*2^n+1$ 的形式。（2^n要大于多项式项数，否则也会炸）

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

解法一

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

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

解法二

考虑拆系数。

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

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

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

于是 $c_1*c_2=(a_1*2^{15}+b_1)(a_2*2^{15}+b_2)=a_1a_2*2^{30}+(a_1b_2+a_2b_1)*2^{15}+b_1b_2$ 。

对这个式子的 $a_1,a_2,b_1,b_2$ 分别做一次FFT，$a_1a_2,a_1b_2,a_2b_1,b_1b_2$ 分别做一次IFFT，总共要用 $8$ 次，常数还是比较大。

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

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

解法三

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

我们现在有四个多项式 $a_1,a_2,b_1,b_2$ ，我们要求 $a_1a_2,a_1b_2,a_2b_1,b_1b_2$ 的值。

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

不妨设 $p=a_1+ia_2$ ， $q=b_1+ib_2$ 。

则 $pq = a_1b_1-a_2b_2+i(a_1b_2+a_2b_1)$ 。

再取 $p'=a_1-ia_2$ ，则 $p'q = a_1b_1+a_2b_2+i(a_1b_2-a_2b_1)$。

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

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

多项式的各种运算

多项式乘法逆（倍增法）

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

我们还是考虑递归。

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

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

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

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

于是我们有：

\[\begin{aligned} f(x)f^{-1}_0(x)&\equiv 1 &\pmod{m'}\\ f(x)f^{-1}(x)&\equiv 1 &\pmod{m'}\\ f^{-1}(x)-f^{-1}_0(x)&\equiv 0 &\pmod{m'} \end{aligned} \]

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

\[f^{-2}(x)-2f^{-1}(x)f^{-1}_0(x)+f^{-2}_0(x)\equiv 0 \pmod{m} \]

两边同乘 $f(x)$ ：

\[f^{-1}(x)-2f_0^{-1}(x)+f(x)f^{-2}_0(x)\equiv 0 \pmod{m} \]

移项得：

\[f^{-1}\left(x\right)\equiv f^{-1}_{0}\left(x\right)\left(2-f\left(x\right)f^{-1}_{0}\left(x\right)\right) \pmod{m} \]

倍增即可。

时间复杂度： $T(n)=T(\frac n2)+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)$，满足以下条件：

$\deg q(x)=n-m$，$\deg r(x)=m-1$
$f(x)=q(x)*g(x)+r(x)$

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

\[f^R(x)=x^{\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)*g(x)+r(x)$ 可知：

\[f(\frac{1}{x})=q(\frac{1}{x})*g(\frac{1}{x})+r(\frac{1}{x}) \]

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

\[f^R(x)=q^R(x)*g^R(x)+x^{n-m+1}r^R(x) \]

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

我们发现上式中 $\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) \pmod {x^{n-m+1}} \]

\[f^R(x)(g^R(x))^{-1}\equiv q^R(x) \pmod {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'=x^{\frac n2}$ 。则有：

\[f(x)-g_0^2(x)\equiv 0 \pmod {m'} \]

两边同时平方，得：

\[\begin{aligned} (f(x)-g_0^2(x))^2&\equiv 0 \pmod m\\ (f(x)+g_0^2(x))^2&\equiv 4f(x)g_0^2(x) \pmod m \end{aligned} \]

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

\[\left(\frac{f(x)+g_0^2(x)}{2g_0(x)}\right)^2\equiv f(x) \pmod m \]

两边同时开方：

\[\begin{aligned} \frac{f(x)+g_0^2(x)}{2g_0(x)}&\equiv g(x) \pmod m\\ 2^{-1}(f(x)g_0^{-1}(x)+g_0(x))&\equiv g(x) \pmod m \end{aligned} \]

倍增计算即可。

时间复杂度： $T(n)=T(\frac n2)+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}}*x^\frac k2$ 即可。

分治FFT

经典例题：给你一个数组 $g_i$ ，求数组 $f_i$ ，使得 $f_i=\sum_{j=1}^if_{i-j}g_j$ 且 $f_0=1$。

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

普通卷积：$h_i=\sum_{j=1}^if_{i-j}g_j$

分治FFT：$f_i=\sum_{j=1}^if_{i-j}g_j$

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

~~但这个性质好像没啥卵用~~

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

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

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

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

\[\begin{aligned} f_3&+=f_2*g_1+f_1*g_2+f_0*g_3\\ f_4&+=f_2*g_2+f_1*g_3+f_0*g_4\\ f_5&+=f_2*g_3+f_1*g_4+f_0*g_5 \end{aligned} \]

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

卷积即可。

时间复杂度： $T(n)=2T(\frac n2)+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}^if_{i-j}g_j \pmod{x^n} \]

由于 $g_0=0$，所以

\[f(x)g(x) \equiv \sum_{i=1}^{n-1} x^i \sum_{j=1}^if_{i-j}g_j \pmod{x^n} \]

并且

\[\begin{aligned} f(x)&\equiv\sum_{i=0}^{n-1} x^if_i\\ &\equiv x^0f_0+\sum_{i=1}^{n-1} x^if_i\\ &\equiv f_0+\sum_{i=0}^{n-1} x^i \sum_{j=1}^if_{i-j}g_j \pmod{x^n} \end{aligned} \]

则有 $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}^\infty \frac{(1-f(x))^i}{i}\\ \]

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

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

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

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

\[\frac 1{f(x)}\cdot f'(x)=g'(x) \]

\[\frac {f'(x)}{f(x)}=g'(x) \]

再同时积分： $$g(x)=\int\left(\frac {f'(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)*f'(x)=g'(x)$

提取系数后就变成了 $g'_n=\sum_{i=0}^{n}g_{n-i}f'_i$

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

即 $ng_n=\sum_{i=1}^{n}(i\cdot g_{n-i}f_i)$，$g_n=\frac 1n\sum_{i=1}^{n}(i\cdot g_{n-i}f_i)$。

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

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

其实，把上述式子再变换也可以求出 $\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\le10^{10^5}$

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

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

\[f(x)^k\equiv f(x)^{k\bmod p}\pmod{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_kx^k+g(x)$，

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

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

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

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

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

快速幂的式子是：

\[f(x)^k\equiv f(x)^{k\bmod (p-1)} \pmod p \]
而在 $n<p$ 的条件下，多项式快速幂的式子是：

\[f(x)^k\bmod {x^n}\equiv f(x)^{k\bmod p}\bmod 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\ne1$ 的情况。

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

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

注意这里要同时保留 $k\bmod p$ 和 $k\bmod (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;
}

多项式牛顿迭代

开个坑，以后再写。

posted @ 2023-02-17 11:45 SqrtSecond 阅读(2895) 评论(0) 收藏举报

刷新页面返回顶部

SqrtSecond

菜死了