多项式:从入门到全家桶
第一次用博客园写学习笔记,写的不好请见谅~
欢迎各位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)\) 上。
定义加法和乘法以后,多项式的基本运算和整数基本一致。(这个大家都会吧)
多项式的导数公式:
多项式的积分公式:
一个长度为 \(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\):
这三个性质在快速傅里叶变换中会被用到。
顺便一提,单位根还有一个比较重要的性质,是: \(\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\) 为例:
我们先按照次数的奇偶来分成两组,然后右边提出来一个 x:
然后再用奇偶次次项数建立新的函数:
那么原来的 \(f(x)\) 就可以表示为:
于是我们可以发现一个性质:它可以分治!
也就是说,我们可以先将原来的式子的奇数项和偶数项分别处理(即分别计算FFT),求出原文中的 \(f_0\) 和 \(f_1\) ,再代入式子求解。
那怎么合并呢?
由于 \(f_1\) 和 \(f_2\) 都只有 \(n/2\) 个点值,那么我们如果用以上那个式子的话,也只能推出 \(f\) 的 \(n/2\) 个点值。
但是,我们可以再考虑以下这个式子:
这样的话,我们能用这个式子求出另外 \(n/2\) 个点值。
你以为这就完了吗?还没有!
观察上面两个式子,我们发现:它只能求互为相反数的 \(n/2\) 对 \(x\) ,而如果在实数范围内考虑,右边就要求 \(n/2\) 个正数的点值,显然行不通。
于是我们考虑单位根。
单位根有很好的的对称性,并且满足平方还是单位根这一特性。(因为 \((\omega_n^k)^2 = \omega_{n/2}^k\) )
如果我们把 \(n\) 个单位根代入,那么式子就会变成这样:
(注意到之前我们让 \(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的过程变成这样:
(就是系数变成它的共轭复数)
然后最终得出来的结果除以 \(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)\) 。
于是我们有:
将最后一个式子平方可得:
两边同乘 \(f(x)\) :
移项得:
倍增即可。
时间复杂度: \(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)\) 即可。
![](https://img2023.cnblogs.com/blog/3051009/202302/3051009-20230221084318342-1008832705.gif)
好吧,真正的多项式除法实际上是带鱼带余除法。
或者说得正式点:给定一个 \(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(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)\) 可知:
两边同时乘上 \(x^n\) 得:
注意到我们只需要求出 \(q(x)\) 和 \(r(x)\) 中的任意一个就能求出另一个。
我们发现上式中 \(\deg q(x)=n-m\) ,而 \(r^R(x)\) 的系数为 \(x^{n-m+1}\) 。
因此我们把式子两边模上 \(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}\) 。则有:
两边同时平方,得:
将 \(4g_0^2(x)\) 移到左边得:
两边同时开方:
倍增计算即可。
时间复杂度: \(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]\) ,则我们需要
于是,\(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);//分治右区间
}
但是先别走!这道题也可以用多项式求逆解决!!!
我们发现
由于 \(g_0=0\),所以
并且
则有 \(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\) 的完整定义在这:
但是,直接用定义式时间复杂度会炸。
首先,对于多项式 \(f(x)\),若 \(\ln(x)\) 存在,则有 \([x_0]f(x)=1\) 。(否则常数项不收敛)
我们知道 \(\ln x\) 的导数为 \(\frac 1x\) ,不妨设 \(\ln f(x)=g(x)\) ,两边同时对 \(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;
}
多项式快速幂
我们考虑封装多项式乘法,再倍增即可。
![](https://img2023.cnblogs.com/blog/3051009/202302/3051009-20230221084318342-1008832705.gif)
好吧,正解是这样的:
我们先保证 \(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}\)
那怎么办呢?还要写个高精?
肯定不可能,我们需要知道一个性质:
这个式子在 \(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;
}
多项式牛顿迭代
开个坑,以后再写。